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Middle  Atmospheric  Periodic-Structure-Associated  Radiance  (MAPSTAR)  is  a 
AFOSR  (2310G201) /Geophysics  Laboratory  sponsored  program.  MAPSTAR  is  a  continuation 
of  a  series  of  ground-based  experiments  designed  to  measure  and  characterize 
time-dependent  mesospheric  emissions  above  large-area  terrains  of  interest  (e.g., 
oceans,  mountains,  and  plains).  Data  from  these  experiments  will  be  used  to 
develop  a  middle-atmospheric  dynamic  structure  data  base  for  Air  Force  modelers  and 
systems  designers.  Space  Dynamics  Laboratory  (SDL)  at  Utah  State  University  (USU) 
participates  in  this  program  under  Geophysics  Laboratory  contract  F19628-87-C-0100. 

The  scientific  report  that  follows  develops  a  model  of  mesospheric  airglow 
observables . 
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CHAPTER  I 


INTRODUCTION 


Introduction 

This  is  a  study  about  the  dynamics  of  the  mesospheric  hydroxyl 
infrareu-emitting  layer.  In  this  chapter,  a  brief  background  of  the 
upper  atmosphere  study  is  presented,  followed  by  a  short  introduction 
to  the  mesospheric  hydroxyl  emission  layer.  Emphasis  then  shifts  to  a 
geometric  consideration  for  the  ground-based  passive  optical  observa¬ 
tion  of  the  layer,  leading  to  a  rarely  studied  effect  in  which 
apparent  oscillations  occur  in  observational  data  which  do  not 
represent  any  physical  property  of  the  layer  but  are  caused  purely  by 
the  geometric  effect.  Finally,  the  objectives  of  this  study  are 
presented  along  with  the  technical  approaches  which  will  be  developed 
in  the  following  chapters. 

Dynamics  of  the  Upper  Atmosphere 

In  the  last  several  decades,  studies  on  chemical  and  thermodynamic 
properties  of  the  upper  atmosphere  have  included  both  physical 
measurements  and  theoretical  modeling.  Our  current  understanding  of 
the  earth's  upper  atmosphere,  derived  from  physical  measurements  to 
date,  is  summarized  in  Fig.  1-1  [1].  The  atmospheric  regions  depicted 
in  this  diagram  are  classified  according  to  altitudinal  variations  in 
temperature.  The  major  focus  of  this  study  is  the  mesospheric  region, 
which  extends  in  altitude  from  60  to  90  km.  In  this  region,  several 
atmospheric  species  exhibit  special  properties  not  found  in  other 
altitude  regions.  For  example,  chemically  reactive 
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species  like  hydroxyl  (OH)  form  an  infrared-emitting  layer  as  a  result 
of  complicated  chemical  processes. 

In  the  first-order  approximation,  the  gravitationally  bound  atmos¬ 
pheric  species  are  in  hydrostatic  equilibrium.  This  approximation 
fails  in  the  stratosphere  and  the  thermosphere  due  to  a  rise  in  local 
temperature,  induced  from  the  generation  of  energetic  species  by  solar 
radiation.  In  the  stratosphere,  the  kinetic  energy  is  produced  by  the 
ultraviolet  photodissociation  of  ozone  according  to  the  reaction  [2], 
13]: 


O3  +  hv  - >  0  +  O2  +  K.E.  (1.1) 

which  provides  the  heat  source.  The  ozone  is  recycled  as: 

0  +  02  +  M  - »  03  +  M  (1.2) 

where  M  is  an  arbitrary  third  body.  In  the  thermosphere,  energetic 
atomic  nitrogen  N  and  atomic  oxygen  0,  together  with  their  ionized 
species,  are  the  major  products.  These  ionized  products  are  more 
chemically  reactive  than  their  neutral  counterparts. 

The  mesopause  is  the  boundary  between  the  turbosphere  and  the 
homosphere,  in  other  words,  the  boundary  between  the  fluid  atmosphere 
below  and  the  molecular  flow  above  [1].  The  mesosphere  is  referred  to 
as  the  region  immediately  below  the  mesopause.  At  the  mesopause, 
abput  86  km  in  altitude,  lies  the  thin  hydroxyl  (OH)  emission  layer, 
which  is  only  a  few  kilometers  thick.  The  formation  of  this  layer  is 
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governed  by  the  photochemistry  and  dynamics  of  active  oxygen  and 
active  hydrogen  species,  including  0,  O3,  H  and  H0£,  as  well  as  by  a 
change  in  pressure  and  temperature  of  the  major  species,  namely  N2 
and  O2.  The  concentrations  and  temperatures  of  the  major  species  can 
be  evaluated  by  the  thermodynamic  equation  of  state,  which  determines 
the  active  species'  concentrations  by  the  cascading  chemical  react¬ 
ions.  These  concentrations  can  be  used  to  calculate  the  amount  of  the 
hydroxyl  with  appropriate  chemical  reactions.  This  process  is 
discussed  in  detail  in  Chapter  II. 

The  Hydroxyl  Emission  Laver 

In  the  night  sky,  the  major  source  of  infrared  emission  is  the 
hydroxyl  (OH)  layer  described  above.  Because  the  lower  atmosphere  is 
transparent  to  portions  of  the  near— infrared  region,  studies  on  the 
upper  atmosphere  have  been  conducted  by  ground-based  optical  observa¬ 
tion  of  infrared  emissions  from  the  layer.  The  fact  that  the  air 
density  in  the  layer  region  is  too  thin  for  aircraft  to  maintain 
altitude  but  too  thick  for  artificial  satellites  to  maintain  orbits 
also  made  the  remote  observation  the  most  common  method  of  studying 
the  dynamics  of  the  mesosphere.  Among  those  mesospheric  studies,  a 
number  of  reports  show  the  existence  of  patchiness  and  wave-like 
structures  in  emissions  from  this  layer  [4],  [5],  (6],  Several 
theoretical  studies  ascertain  the  formation  mechanisms  of  these 
spatial  structures  [7].  Evidence  has  accumulated  in  support  of  the 
theory  that  the  so-called  atmospheric  internal  gravity  wave  is  the 

l 

principal  wave -formation  mechanism  in  many  cases  [5],  [6].  A  detailed 
theory  of  the  atmospheric  gravity  wave  is  presented  in  Appendix  A. 
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The  consensus  now  is  that  the  internal  gravity  wave  is  also  important 
in  the  dynamics  and  heating  of  other  mesospheric  species,  including 
the  turbulent  transport  of  0  and  O3  [8]. 

Wave  Structures  in  the  Mesospheric 
Hvdroxvl  Emission  Laver 

An  observable  from  the  mesospheric  hydroxyl  layer  is  the  spectral 
radiance  of  the  infrared  emissions  from  which  the  rotational  tempera¬ 
ture  can  be  calculated.  A  detailed  theory  of  the  rotational  tempera¬ 
ture  is  given  in  Appendix  B  [  9  J .  The  OH  rotational  temperature  is 
assumed  to  be  representative  of  the  kinetic  temperature  in  the  meso¬ 
sphere  [10],  Both  the  calculated  rotational  temperature  and  the 
emission  intensity  of  the  OH  layer  show  temporal  variations.  These 
variations  often  exhibit  ordered,  wave-like  forms  which  are  functions 
of  time  and  space. 

In  early  studies  of  mesospheric  hydroxyl  dynamics,  a  parameter 
usually  denoted  as  17  (eta)  was  introduced  by  Krassovsky  [4]  and  used 
by  others  [7],  [11]  to  characterize  the  sinusoid  oscillations.  The  rj 
is  defined  as: 


*} 


_  dl/I 
dT/T 


(unitless) 


(1.3) 


where  T  is  the  stationary  temperature,  I  is  the  stationary  emission 
intensity  and  dT  and  dl  are  their  oscillation  amplitudes.  The 
absolute  value  of  17,  first  evaluated  to  have  a  value  of  about  5  j4j, 
has  recently  been  extended  to  a  complex  number  based  upon  model 
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calculations  (7],  [11].  This  extension  allows  phase  shifts  between 
the  oscillations  i-  temperature  and  emission  intensity  to  be  taken 
into  account.  This  phase  shift  has  been  observed  in  some  cases  [12]. 

In  this  paper  presents  a  possible  mechanism  for  the  phase  rela¬ 
tionship  between  the  temperature  and  intensity.  The  first  approach  is 
to  show  that  the  geometrical  observational  set-up  provides  a  relative 
phase  shift  between  intensity  and  temperature  oscillations.  Under 
typical  conditions,  the  temperature  oscillation  contains  a  spurious 
frequency  component,  an  observable  oscillation  which  is  an  artifact  of 
the  geometric  effect.  The  second  approach  uses  a  set  of  chemical  and 
thermodynamical  equations  to  model  the  layer.  The  theory  of  the 
internal  gravity  wave  is  included  to  describe  the  characteristics  of 
the  wave  itself  as  well  as  its  effect  upon  the  OH  layer.  This 
approach  is  based  upon  those  used  in  several  previous  studies  [11], 
[13].,  The  two  methods  can  be  combined  and  extended  into  more  general 
modeling,  allowing  simulations  of  observed  emission  oscillations  with 
off-zenith  observational  configurations.  This  model  extension  is 
unique  to  this  study  and  will  serve  to  facilitate  many  practical 
experimental  setups.  Based  upon  such  model  evaluations,  recommenda¬ 
tions  for  experimental  setups  which  will  achieve  better  scientific 
understanding  of  the  mesospheric  dynamics  will  be  presented. 

It  is  assumed  in  the  definition  of  rj  that  the  wave  is  monochro¬ 
matic  sinusoid.  In  Chapter  III,  such  assumption  is  tested.  It  is 
found  that  the  dynamic  model  predicts  that  the  wave  frequency  varies 
as,  a  function  of  time.  This  finding  is  used  to  construct  the  signal 
processing  system  in  Chapter  V. 
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The  concept  of  q  was  Introduced  using  simple  data  of  low  temporal 
resolution.  It  Is  not  obvious  which  quantity  should  be  used  to 
calculate  *7  in  the  higher  resolution  data  available  today.  Several 
mathematical  considerations  and  tools  are  introduced  to  process  such 
data  in  a  way  useful  for  theoretical  model  evaluation.  Some  discuss¬ 
ion  about  q  is  presented  in  Chapter  VI . 

Scope  and  Objectives 

The  specific  objectives  of  this  study  are  summarized  as  follows: 

1.  To  discuss  Krassovsky's  theoretical  17  using  linearized  chemical 
equations  and  dynamic  internal  gravity  wave  models,  taking  into 
account  the  geometric  effect. 

2.  To  develop  a  procedure  to  simulate  off-zenith  observations  using 
the  model  for  zenith  observations  with  geometric  effect  modifica¬ 
tions. 

3.  To  relate  the  tj  calculated  from  theoretical  models  with  the  values 
calculated  from  the  measurement  data,  taking  into  account  both 
magnitude  and  phase  of  oscillations. 

4.  To  formulate  an  appropriate  signal  processing  technique  that  makes 
possible  a  reasonable  estimation  of  gravity  wave  parameters  from 
observational  data. 

5.  To  perform  an  evaluation  of  dynamic  hydroxyl  emission  layer 
models . 

6.  To  determine  the  relevancy  of  the  concept  of  a  complex  q  and  the 
>  associated  dynamic  models. 

7.  To  consider  the  possibility  of  improving  the  acquisition  of 
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knowledge  of  the  hydroxyl  layer  dynamics  and,  if  such  improvement 
is  found,  to  make  appropriate  recommendations  for  further  study. 
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CHAPTER  II 

BASIC  MODEL  FORMULATION 


This  chapter  describes  a  detailed  development  of  the  chemical 
model  and  a  brief  introduction  of  the  overall  picture  of  the  problem. 

Emission  Layer  Formulation 

In  this  section,  a  set  of  chemical  and  photochemical  reactions 
based  upon  an  equation  set  with  the  updated  reaction  rate  constants 
used  for  ozone  model  calculations  is  presented  [3],  [14].  As 

introduced  in  Chapter  I,  the  mesopause  is  the  narrow,  upper  region  of 
the  mesosphere  in  which  the  ratio  of  the  major  odd  oxygen  consti¬ 
tuents,  atomic  oxygen  and  ozone,  changes  [3],  [15],  The  atomic  oxygen 
is  created  by  photo  dissociation 

02  +  hi/  - >  0  +  0  (2.1) 

during  the  daytime.  The  ozone  is  created  by  the  reaction  [3]: 

02  +  0  +  M  - »  O3  +  M  (2.2) 

where  M  is  a  third  body  such  as  N2,  02  or  0.  The  atomic  oxygen 
created  by  reaction  (2.1)  is  the  major  constituent  of  the  upper 
atmosphere  above  220  km  [1],  Table  I  summarizes  typical  values  of 
physical  importance  in  the  mesospheric  region  [1],  [3].  As  can  be 
evaluated  from  the  values  given  in  Table  I  and  using  equation  (2.2), 
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TABLE  I 

Typical  Physical  Quantities  in  Mesosphere  [3], 
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Species  Concentration  (cm  3) 


n2 

10U 

02 

1013 

0 

1010 

o3 

108 

Temperature  150  K  (Summer)  220  K  (Winter) 


the  atomic  oxygen  has  a  lifetime  of  10  to  15  hours  at  85  to  95  km. 
Consequently,  the  concentration  of  this  species  can  be  regarded  as 
constant,  even  during  the  nighttime. 

The  ozone-creation  profile  resulting  from  (2.2),  which  is 
identical  to  its  concentration  profile,  has  a  sharp  peak  because  the 
creation  mechanism  of  ozone  is  similar  to  the  formation  of  the 
Chapman  layer  which  has  a  gaussian-type  distribution.  This  ozone  peak 
occurs  in  the  stratosphere,  at  an  altitude  of  about  50  to  70  km  [1]. 


.  —  /O  ON 

-/111  \  L.  .  4.  )  Cllfc 


;nc  occurs  with  the  destruction  of  atomic 


oxygen.  Therefore,  there  is  a  transition  of  the  major  odd  oxygen 
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component  from  ozone  to  atomic  oxygen  with  increasing  altitude.  This 
transition  occurs  in  a  narrow  region  at  the  mesopause  [3].  In  this 
region,  one  particular  reaction  is  the  major  contributor  to  the 

formation  of  the  hydroxyl  layer  [3].  (15).  This  reaction  is  the 

hydration  of  ozone  [15] 

H  +  03  - *  OH*  +  O2  (2.3) 

The  vibrationally-excited  hydroxyl  in  the  ground  electronic  state, 
denoted  as  OH*,  can  emit  photons.  The  wavelengths  of  the  resultant 
photon  emissions  fall  in  the  infrared  and  visible  spectral  regions. 
Since  the  hydroxyl  concentration  is  determined  by  the  balance  of  the 
atomic  oxygen  and  the  ozone,  both  of  which  are  highly  chemically 
active,  the  hydroxyl  forms  into  a  relatively  thin  layer  at  the 

transient  region  of  the  odd  oxygen  species.  Most  investigators  report 
the  infrared  emission  region  to  be  a  relatively  thin  layer,  usually 
about  7  km  thick  [10],  consistent  with  the  dynamic  model. 

Formation  of  the  hydroxyl  airglow  layer  can  be  numerically 
simulated  using  a  computer;  that  is,  successive  values  of  the 
concentrations  of  the  species  of  interest  can  be  calculated  by 
employing  a  small  time  step  dt  in  the  equations.  This  approach  is 
explained  in  detai?  in  Chapter  III.  Simplifications  of  the  process 
are  discussed  in  the  following  section.  Table  II  summarizes  the  most 
important  chemical  reactions  related  to  the  mesospheric  hydroxyl 
airglow  layer  [3], 
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TABLE  II 

Chemical  Reactions  Related  to  the  Formation  of  Hydroxyl  Layer  [3). 
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REACTION 

NUMBER  as  in  (3) 


REACTION 
COEFFICIENT 
T:  temperature 

(K) 


R13  0  +  02  +  0 
R14  0+02+02 
R15  0  +  02  +  N2 
R17  0  +  OH 
R29  H  +  O3 


— >  O3  +  0 
— >  O3  +  O2 
— >  O3  +  N2 
— >  02  +  H 
— >  OH*  +  02 


2.15E-34  e(345/T) 
2.15E-34  e(345/T) 
8.82E-35  e(575/T> 
2.3E-11  e(*9°/T) 
1.4E-10  e(-27°/T> 


Chemical  Model 
Simplified  Model 

In  the  simplified  chemical  model,  the  principal  reactions  are  as 
follows  [3] : 


H  +  03 
OH*  +  M 

OH* 


Si 


->  OHJ  +  02 


'1J 


>  OHJ  +  M 


■ij 


■>  OHj  +  hi/£j 


(2.4) 

(2.5) 

(2.6) 


where  i . j  -  0...9  are  vibrational  states  and  0H|  represents  the 
hydroxyl  which  is  vibrationally  excited  to  state  i.  Equation  (2.5) 
represents  the  quenching,  a  transition  between  states  without  photon 
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emission.  Equation  (2.6)  represents  the  photo  emissive  transitions 
which  give  rise  to  the  mesospheric  airglow  layer.  The  coefficients 
ljj  are  known  as  Einstein  coefficients  [10].  The  following  set  of 
equations  follow: 


d[OH*] 

dt 

-  gi[H]{03] 

(2.7) 

d  ( OH* ] 
dt 

i-1 

-  -  ^2^  k£j [M] [OH*] 

(2.8) 

d[OH*] 

dt 

i-1 

2  lij  [OH*] 
j-0 

(2.9) 

This  set  can  be  represented  by  a  10  by  11  matrix;  however,  with 
the  following  assumptions,  most  of  the  off-diagonal  elements  can  be 
neglected: 

1]  Production  of  OH*  occurs  only  for  i-7,8  and  9. 

2]  Quenching  shown  by  (2.5)  occurs  only  when  i-j  -  1. 

These  assumptions  are  represented  by  the  following  relationships: 

gi  H  0  (i>6) 

-  0  (otherwise)  (2.10) 


kij  ^  0  (i-j-1) 

-  0  (otherwise)  (2.11) 

With  these  constraints,  equations  (2.7),  (2.8)  and  (2.9)  are 

reduced  to  the  following  set: 
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d[OHf] 

dt 

i-1 

-  -  [OH*]  (  2  (lij)  +  kii-iCM]  )  +  gi[H]  [03] 

j-o 

(2.12) 

This  matrix  equation  can  be  expressed  as  follows: 

& 

dt 

-  AX  +  B 

(2.13) 

where 

i-1 

A  -  (  2  (In)  +  kii.x  [M]  ) 

j-o 

B  -  (  gi(H] [03]  ) 

X  -  (  [OH*]  ) 

(2.14) 

This  has  a  solution  of: 

X  -  exp[-/A(t)dt  +  exp(— jA(t)dt)  ] 

./exp(/A(t)dt)  B(t)dt  X(0)  (2.15) 

where  X(0)  is  the  temporrl  initial  values  set.  Equation  (2.15)  can  be 
reduced  to  the  following  form  with  (2.13): 

X  -  exp  (At)  X(0)  -  exp  (At)  A_1B  +  A*XB  (2.16) 

where  A“^  represents  the  inverse  matrix  of  A. 

Further  simplification  is  possible.  Previous  studies  have  shown 
that  the  reactions  (2.8)  and  (2.9)  have  time  constants  on  the  order  of 
tens  of  milliseconds.  The  observational  data  are  discrete  time  series 
with  tens  of  seconds  time  steps.  Because  this  is  almost  one  thousand 
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times  longer  than  the  time  constants  of  (2.8)  and  (2.9),  the  con¬ 
centration  of  excited  hydroxyl  of  all  levels,  [OH*]  i-1-9,  can  be 
approximated  to  be  at  equilibrium  for  all  data  points.  This  is 
equivalent  to  taking  the  exp(At)  in  (2.16)  as  zero  which  gives  the 
following  simple  solution  for  X: 

X  -  A_1B  (2.17) 


Elements  of  A  and  B 

The  elements  of  the  transition  matrix  A,  defined  in  the  previous 
section,  are  deduced  from  the  Einstein  coefficients  for  the  excited 
hydroxyl  (the  emissive  transition  rates  among  the  excited  states)  and 
the  quenching  rates  for  each  of  the  states.  These  numbers  are 
provided  in  Table  III  [15].  Note  that  the  matrix  A  is  upper-half 
triangular  with  non-zero  diagonal  elements. 

The  elements  of  the  column  vector  B  are  the  reaction  rates  g^  of 
equation  (2.4),  which  have  non-zero  values  only  for  i-7,8,9.  However, 
a  recent  study  indicates  that  the  elements  of  B  may  be  non-zero  for 
other  levels  as  well,  including  i-4-9  [15].  These  values  are  also 
included  in  Table  III. 

Model  Extension 

The  simple  model  described  in  the  previous  section  is  extendable 
according  to  the  following  rules,  allowing  for  the  formulation  of 
appropriate  models  under  different  assumptions-:  first,  to  add  an 
additional  chemical  reaction  to  the  model,  there  are  three  cases  which 
depend  on  whether  the  OH  itself  is  involved  in  the  reaction. 
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TABLE  III 

Constants  Used  in  the  Model 


16 


Quenching  Rate  *10‘ 1 3 

v  -  1  2  3  4567  8  9 

0.000  0.003  0.23  1.1  3.5  6.5  9.6  11.5  8.9 

Generation  Ratio  (%) 

v  -  1  2  3  4  5  6  7  8  9 


0.0  0.0  0.0  0.0  0.0  0.0  15.0  31.0  38.0 


Einstein  Coefficient 


V' 

1 

2 

3 

4 

5 

6 

7 

8 

9 

#  — 

0 

-31.7 

40.6 

67.9 

13.5 

2.81 

0.59 

0.13 

0.03 

0.01 

1 

0.0 

-67.0 

35.0 

116.3 

28.5 

7.39 

1.82 

0.46 

0.13 

2 

0.0 

0.0 

-106.8 

22.1 

165.9 

47.4 

14.9 

4.43 

1.24 

3 

0.0 

0.0 

0.0 

-152.4 

8.33 

213.5 

68.5 

24.8 

8.90 

4 

0.0 

0.0 

0.0 

0.0 

-205.6 

0.31 

256.4 

90.7 

36.9 

5 

0.0 

0.0 

0.0 

0.0 

0.0  - 

•269.2 

5.12 

291.8 

110.2 

6 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0  ■ 

■346.9 

31.1 

322.8 

7 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

-443.3 

90.5 

8 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

-570.6 

Case  1:  OH*  is  in  the  product  only. 

hi 

X  +  Y  - >  OH*  +  Z  (2.18) 


where  i  is  the  appropriate  vibrational  level.  This  equation  leads  to 
the  following  equation: 
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d[OH*] 

— ^  IX]  [Y]  (2.19) 

The  modification  is  performed  on  the  vector  B  by  adding  the  term 
hi[X][Y]  to  the  i-th  row. 

hi  -bi+hi(X][Y)  (2.20) 


Case  2:  OH*  is  both  in  the  reactant  and  in  the  product. 


X  +  OH* 


OHj  +  Z 


(2.21) 


which  leads  to  the  following  equation  set: 


d[OH*] 

dt 


-  hij  [X] [OH*] 


d[OH* ] 
dt 


hy  [X]  [OH*] 


(2.22) 


(2.23) 


The  modification  is  then  performed  on  the  matrix  A  by  adding  — hi j [ X ] 
to  its  ajj  element: 

a;j  -  a{  j  -  h| j  [X]  (2.24) 


Case  3:  OH*  is  in  the  reactant  only. 
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Hj 

X  +  OHj  - >  Y  +  Z  (2.25) 

which  leads  to  the  following  equation: 
d[°H*f] 

— -  ~hj  [X]  [OHj  ]  (2.26) 


This  is  accomplished  by  adding  — h j  [ X )  to  the  ajj  element  of  the  macrix 
A: 


a. .  -  a. .  -  h: [X]  (2.27) 

jj  n  J 

Using  these  techniques,  the  previous  model  can  be  simplified  as 
follows.  This  simpler  model  is  appropriate  when  the  time  scale  of 
interest  is  on  the  order  of  hours,  but  may  suffer  some  problems  when 
the  time  scale  is  less  than  30  minutes.  The  reaction  R29  has  a  time 
constant  of  about  10  minutes.  Therefore,  when  the  time  scale  of 
interest  is  much  larger  than  lO  minutes,  on  the  order  of  hours,  the 
concentration  of  ozone  can  be  regarded  at  its  equilibrium  value.  In 
such  a  case,  the  reactions  R13,  14  and  15  can  be  substituted  as  the 
production  of  OH  instead  of  R29  [15].  This  simplifies  the  model 
because  the  reliable  observational  data  on  [O3]  and  [H]  are  not 
readily  available.  In  this  case,  vector  B  in  the  equation  (2.14)  is 
replaced  by: 

B  -  (  gi [0] [02 ] [M]  )  (2.28) 
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where  g'  denotes  the  appropriate  coefficient  and  [M]  is  either  02,  0 
or  N2,  depending  upon  which  reactions  are  to  be  used.  The  complete 
expression  is  provided  in  Table  IV. 


Thermodynamic  Equation 


The  thermodynamic  equation  of  state  is  represented  as: 


(2.29) 


where  N  is  the  number  of  molecules  in  volume  V,  p  is  the  pressure,  R 
is  the  gas  constant,  8.314  J/mol  K,  and  T  is  the  temperature.  Under 
the  adiabatic  condition,  this  equation  reduces  to  [4]: 


(2.30) 


where  7  is  the  ratio  of  the  constant-pressure  specific  heat  Cp  to  the 
constant- volume  specific  heat  Cy  and  K  is  a  constant.  For  future 
uses,  the  following  formula  provides  an  alternative  expression  for  the 
adiabatic  state  equation: 

p-r-1  -  VXL  (2.31) 

m 


where  k'  is  a  constant  of  proportionality,  m  is  the  average  molecular 
mass  of  air  and  p  is  the  average  mass  density  of  the  air  (kg/cm3). 
The  expression  (2.31)  is  useful  for  detailed  modeling  because  it 
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TABLE  IV 

Matrix  A  and  Vector 
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MO  0  S4 10|  |Oj  1(H)  f,IOUO,J(H|  tolOltOjUM] 


21 


involves  only  Che  observable  quantities  in  the  mesopause  region. 
Internal  Gravity  Wave 

Evidence  suggests  that  the  so-called  atmospheric  internal  gravity 
wave  (IGW)  is  the  source  of  spatial  and  temporal  oscillations  observed 
in  hydroxyl  infrared  emissions  in  the  mesospheric  layer  [5],  [16], 

[17].  Hines  summarized  linearized  IGW  model  (Appendix  A)  [13],  [18]. 
Major  characteristics  of  IGW  related  to  study  of  the  mesospheric  layer 
can  be  summarized  as  follows: 

1.  IGWs  are  dispersive.  This  means  a  wave's  packet  and  phase 
velocity  depend  upon  its  oscillation  frequency.  This  property 
sets  an  upper  limit  to  the  oscillation  frequency  that  can  reach 
to  mesospheric  altitudes. 

2.  The  propagation  speed  of  the  IGW  depends  upon  its  propagation 
direction  and  the  frequency,  suggesting  that  the  observed  IGW  at 
a  certain  height  is  almost  monochromatic  [13],  [18]. 

3.  The  temperature  inversion  below  the  mesosphere  prevents  short- 
period  components  of  IGW  originating  at  lower  altitudes  from 
propagating  upward  further  into  the  mesosphere.  This  makes  the 
upper  limit  of  allowed  frequencies  of  IGWs  in  the  mesosphere 
smaller  than  the  Brunt -Vaisala  frequency. 

4.  Energy  loss  due  to  viscosity  becomes  significant  for  IGWs  with 
short  wavelengths . 

5.  The  phase  velocity  and  group  velocity  are  .almost  perpendicular 
to  each  other.  Therefore,  an  IGW  originated  at  a  lower  altitude 
and  propagating  upward  to  the  mesopause  region  at  group 
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velocity,  appears  locally  to  propagate  downward  at  phase 
velocity.  Most  observations  measuring  only  a  few  periods  of 
oscillation  are  likely  to  observe  the  phase  velocity,  or  "phase 
trace  speed"  as  Hines  defined  [18]. 

These  properties  of  IGWs  put  limits  on  the  characteristics  of  a 
particular  IGW  that  can  reach  in  the  mesospheric  OH  airglow  layer 
region,  such  as  the  oscillation  period,  horizontal  and  vertical 
wavelengths  and  propagation  direction.  Table  V  provides  equations 
regarding  properties  of  an  IGW  in  general  as  well  as  additional  equa¬ 
tions  that  model  the  characteristics  of  IGWs  in  the  mesospheric  region 
[13],  [18],  [19].  Fig.  2-1  depicts  the  general  propagation  property 
of  IGWs  described  by  those  equations  [13].  The  family  of  ellipsoids 
represents  the  acoustic-wave  solutions  to  the  equations  in  Table  V. 
The  family  of  hyperbolas  represents  the  IGWs.  Each  curve  is  marked  by 
oscillation  periods  in  minutes.  In  this  example,  7  is  1.4,  g  is  9.5 
(m/s2),  uA  corresponds  to  4.4  minutes  and  wg  is  4.9  ">inutes.  The 
expected  characteristics  of  an  IGW  in  tha  mesospheric  region  are  shown 
in  Fig.  2-2  [13],  [18].  Fig.  2-2(a)  and  (b)  depict  the  parameters  of 
the  IGW's  that  can  propagate  through  the  atmosphere  between  the  tropo¬ 
sphere  and  mesosphere  [13].  Fig.  2-2(c)  depicts  the  i  'ationship 
between  the  oscillation  period,  group  velocity,  propagation  elevation 
angle  and  phase  velocity  of  an  IGW  [18],  Fig.  2-2(a)  illustrates  the 
spectrum  of  wavelengths  available  at  90  km,  following  propagation  from 
the  lower  atmosphere,  shown  by  the  unhatched  area  [13].  Fig.  2-2 (b) 
depicts  propagation  modes  at  meteor  heights.  The  periods,  measured  in 
minutes,  are  shown  in  boxes  on  the  corresponding  constant-period 
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-1.  Propagation  property  of  IGW  [13J 


TABLE  V 

Dynamic  IGW  Equations  [13] 


J*  -  w2c2(k2+k2)  +  (7— 1 ) §2k2  - 


2  2  2 
7  g  w 
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TABLE  V 

Dynamic  IGW  Equations  [13]  (continued) 


wA  -  7g/2c 

wg  -  Vy-l  g/c 


< , ?■  -  c^(k^  +  k|)  Acoustic  Wave 


Fig.  2-2.  Propagation  modes  of  IGW  [13],  [18] 
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Fig.  2-2.  Propagation  modes  of  IGW  [13],  [18]  (continued) 


contours  (solid  lines).  The  limits  of  the  permitted  spectrum,  as 
determined  by  viscous  damping,  are  shown  for  several  heights  by 
broken  curves.  Modes  lying  above  and  to  the  right  of  these  curves 
are  excluded.  The  modes  subject  to  reflection  are  shown  by  dotted 
lines.  Modes  lying  below  these  curves  cannot  proceed  from  the  lower 
atmosphere  to  the  upper  [13].  Fig.  2-2(c)  depicts  an  IGW  front 
(bold  curve)  in  the  course  of  expansion  from  a  point  source  at  the 
origin  of  coordinates  [18].  It  is  depicted  at  unit  time  after  its 
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emission,  if  the  spatial  coordinates  are  taken  to  be  measured  in 
units  such  that  the  speed  of  sound  is  1.  The  wave  period  at  the 
front  is  indicated  by  multiples  of  Tg«2?r/wg  at  various  points  along 
the  front.  Inclined  line  segments  behind  the  advancing  front 
indicate  planes  of  constant  phase  immediately  behind  the  front.  The 
bulk  of  some  observable  disturbances  may  proceed  more  slowly  than 
the  frontal  speeds  and  may  be  confined  at  this  particular  moment  of 
time  to  a  region  such  as  the  shaded  region  [18].  As  a  first-order 
approximation  that  the  hydroxyl  layer  is  horizontal  with  a  constant 
height,  the  pattern  can  be  represented  by  the  intersection  of  this 
curve  and  the  layer  plane.  In  a  more  accurate  approximation,  one 
must  take  into  account  the  chemical  reactions  governing  the  hydroxyl 
concentration  involving  the  odd  oxygen  constituents  whose  concentra¬ 
tions  are  first  changed  by  the  incoming  IGW. 

Geometric  Effect 

Because  this  study  involves  a  comparison  of  theoretical  models 
with  ground-based  optical  observational  data,  the  geometric  effect 
must  be  taken  into  account.  The  problems  associated  with  observa- 
tionalgeometry  has  been  pointed  out  briefly  [11]  but  not  worked  out 
in  great  detail.  The  approach  will  be  illustrated  here  using  a 
simplified  model.  A  more  elaborate  model  will  be  introduced  in 
chapter  IV. 

In  ground-based  optical  measurements,  the  instrument  is  often 
oriented  to  look  up  at  the  sky  in  the  near  zenith,  resulting  in  a 
geometric  configuration  like  the  one  shown  in  Fig.  2-3.  Since  the 
OH  airglow  layer  is  about  5  to  10  km  thick,  the  observation  is  not 
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Fig.  2*3.  Ground-based  observation  geometry 

an  emission  from  a  localized  small  volume  but  rather  of  the 
integrated  volume  emission  within  the  field  of  view  along  the  line 
of  sight.  The  combination  of  a  linearized  IGW  model  (13]  «vid 
experimental  airglow  data  [4]  can  be  used  to  demonstrate  the  sig¬ 
nificance  of  this  geometrical  effect. 

A  simple  linearized-IGW  theory  includes  a  directional  disper¬ 
sion  relation  (20): 

o  ”  cos  6  (2.32) 


where  w  is  the  frequency  and  6  is  the  angle  of'  the  packet  velocity 
relative  to  the  zenith  and  n  is  the  Brunt-Vaisala  frequency  [20]. 
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The  phase  velocity  is  perpendicular  to  the  packet  velocity,  that  is, 
perpendicular  to  the  phase  plane.  Table  VI  and  Fig.  2-4  demonstrate 
several  values  of  6  with  typical  values  of  w.  Table  VII  illustrates 
several  characteristics  for  an  airglow  observation  (5).  This 
example  shows  an  IGW  with  a  period  of  11  minutes  and  a  wavelength  of 
14  km  in  the  vertical  direction.  If  the  wave  is  represented  by: 

I(t,z)  -  I0  +  AI  sin(ut-kj)  (2.33) 

where  I0  is  a  constant,  or  a  DC  component,  and  AI  is  the  amplitude 
of  the  oscillation,  then  the  observed  airglow  intensity  is: 

Ah 

lobs ( t)  -  /  I(t,z)dz 
0 

Ah 

-  I0Ah  +  AI  /  sin(wt-kz z)dz  (2.34) 

0 

instead  of: 

I  ,  (t)  -  I0Ah  +  AIAh  sin(wt-k, z)  (2.35) 

ref  * 


TABLE  VI 

Wave  Period  and  Slant  Angle  6 


Period  (min) 

10 

20 

30 

60 


6  (degree) 

60.0 

75.5 

80.4 

85.2 
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TABLE  VII 

Observed  Wave  Parameters  [5] 


30 


Ax  26.0  ±  0.3  km 

vx  38.0+2.0  m/s 

Az  14.5  ±  2.0  km 

vz  2 1 . 0  ±  3 . 0  m/s 

—  11 .4  ±  0.6  min 


From  (2.32),  d  is  about  57  degrees  for  a  typical  value  of  0, 
the  Brunt-Vaisala  frequency.  If  the  layer  is  approximated  by  a  5km- 
thick  horizontal  slab,  the  vertical  line  of  sight  will  integrate 
about  5/14  of  one  oscillation  period,  or  about  128  degrees  of  a  360 
degrees  phase.  Under  these  conditions,  the  integration  of  (2.34) 
is : 


I  (t)  -  I0 Ah  +  ^  (2(l-cos(kzAh))]1/Z  sin(wt+>&)  (2.36) 

obs  Kz 

*  -  arcCan(£osIkiAh}zi)  .  _  kl4h  (2 .37) 

sin(kzAh)  2 


Fig.  2-5  is  a  plot  of  the  apparent  amplitude  and  phase  as 

functions  of  the  thickness  of  the  emitting  layer.  Ac  Ah-5  km,  is 

about  -55  degrees  and  the  amplitude  is  about  85%  of  I  .  Since  Ah 

ref 

is  usually  unknown  and  therefore  needs  to.be  estimated,  the 
geometric  effect  can  cause  considerable  uncertainty  in  both  the 
amplitude  and  phase  of  the  apparent  emissions.  The  effect  could 
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mask  out  amplitude  and  phase  modifications  caused  by  other  factors 
such  as  chemical  reaction- rate  dependencies.  Therefore,  further 
study  of  the  geometric  effect  is  needed  in  order  to  model  mesosphe¬ 
ric  dynamics  and  to  compare  the  model  with  experimental  data. 

Similarly,  it  is  possible  to  assess  the  geometric  effect  on  the 
apparent  temperature  of  the  layer.  The  rotational  temperature, 
calculated  from  observational  data,  is  assumed  to  be  equal  to  the 
local  kinetic  temperature.  The  actual  observed  value  is  "bright¬ 
ness-weighted  temperature"  through  the  layer  along  the  line  of 
sight.  In  each  of  three  cases  depicted  in  Fig.  2-6,  e  tempera¬ 
ture  is  assumed  to  be  a  constant  at  a  same  altitude.  The  tempera¬ 
ture  at  height  h  is  approximated  by  the  following  equation: 


T(h)  -  T 


«-  T i 


Ah 


<h-h0) 


(2.38) 


The  observed  temperature,  the  brightness -weighted  average  tempera¬ 
ture,  is  therefore: 


Tobs ( t) 


ho+Ah 

J"  T(h)I (h)dh 

ho 


ho  +Ah 

/  J  1 (h)dh 

ho 


(2.39) 


This  integral  can  be  evaluated  as: 


Tobs ( t) 


!o Ah  ■CXMX+Tmtn)  +  Asin(fc>t+»)' 
IoAh  +  Bsin(wt+<l>) 


(2.40) 
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where 


A2 


2  2 
T  +  T  . 

max  min 


+ 


2  AT2 
(kz  Ah)2 


(l-c^stkz Ah) ) 


+ 


2AT2 
kj  Ah 


sin(kzAh)  -  2T 

max 


T  . 

mi  n 


cos(kz Ah) 


(2.41) 


-1  Tmaxkz  Ahcos (kz  Ah)— Tmi nkz Ah— ATsinfkz Ah) 
an  Tmaxkz  Ahsin(kz Ah)+AT(cos (kz Ah)-1) 


(2.42) 


AT 


T  -T  . 
max  mi  n 


Fig.  2-7  shows  a  family  of  curves  of  the  function: 


(2.43) 


Y  -  ,/■— ±— S-lsIIX.  (2.44) 

5  +  sin(T+$) 

where  4  is  taken  as  0,  30,  60,  90  and  180  def~  »«■  Since  (2.44)  is 

a  generalized  normal  form  of  (2.40),  this  illustrates  the 

form  of  T  (t)  for  several  phase  shifts  betw>_  v  and  4.  Note  that 

in  the  presence  of  a  sinusoidal  IGW,  even  when  there  is  no  local 

oscillation  in  temperature,  the  observed  temperature  shows  non- 

sinusoidal  oscillation  at  the  same  frequency  as  the  IGW.  Case  2  is 

the  same  equation  with  T  and  T  .  interchanged.  Case  i  can  be 

max  min 

considered  a  combination  of  the  two  cases  as  shown  in  Fig.  2-8. 
The  importance  of  the  geometric  effect  is  that  even  without  any 
local  temperature  fluctuations,  the  airglow  OH  rotational  tempera¬ 
ture  shows  an  apparent  oscillation  due  to  the  temperature  difference 
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Fig.  2-8.  Relation  between  temperature  models 

with  respect  to  altitude  within  the  emitting  layer.  This  suggests 
that  when  there  is  a  local  temperature  change  due  to  an  IGW,  such  a 
change  will  be  observed  with  the  geometric  effect  "modulation." 
Therefore,  it  is  necessary  to  consider  the  geometric  effect  when 
analyzing  temperature  distributions  with  respect  to  altitude. 

The  previous  simple  examples  present  the  general  geometric 
effect;  however,  several  assumptions  were  made  in  Fig.  2-6  which 
could  be  too  limiting.  The  first  such  assumption  is  that  the 
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emission  layer  is  approximated  by  an  uniform  horizontal  slab.  This 
is  equivalent  to  assuming  that  the  layer  has  a  uniform  volume 
emission  rate  throughout  an  entire  layer  with  sharp  top  and  bot  om 
boundaries . 

The  second  assumption  is  that  the  layer  is  time  invariant  while 
the  IGW  passes  through.  This  is  equivalent  to  ignoring  both 
temperature  and  concentration  changes  due  to  an  adiabatic  process. 
Because  some  chemical  models  of  the  layer  are  available,  including 
the  model  using  MSIS  [21],  more  detailed  simulation  of  the  geometric 
effect  will  be  presented  in  chapter  IV. 

Krassovsky's  Parami  • .  ~  Eta 

This  section  discusses  the  parameter  introduced  by  Krassovsky 
[4]  and  since  used  in  several  mesospheric  infrared  emission  studies 
[7],  [11],  denoted  with  the  character  r\  (eta)  of  the  Greek  alphabet. 
An  attempt  will  be  made  to  extend  the  concept  of  eta. 

1.  First  Definition  of  eta 

Krassovsky's  definition  of  eta  is  (4): 


AI /I 
AT/T 


(unitless) 


(2.45) 


where  I  is  the  volume  emission  rate,  T  is  the  temperature  and  AI  and 
AT  are  the  respective  oscillation  amplitudes.  Theoretical  works  to 
explain  the  value  and  the  temporal  behavior  of  eta  with  regard  to 
the  physical  conditions  of  the  hydroxyl  layer'  can  be  generalized 
into  two  major  groups.  The  first  group  is  represented  by  the  work 
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of  Krassovsky  [4] .  In  this  model  it  is  assumed  that  the  adiabatic 
change  of  an  air  parcel  is  caused  by  internal  gravity  waves.  Here, 
the  vibrationally  excited  molecule  X02  is  assumed  to  be  formed 
according  to  the  reaction: 

X  +  02  +  M  - >  X02  +  M  (2.46) 

where  X  denotes  either  0  or  H  atoms.  The  reaction  rate  coefficient 
Rx  is  expressed  as: 

Rx  -  a  T“n  (2.47) 

where  the  value  of  the  exponent  n  depends  on  each  particular  X. 
Krassovsky  gave  n-4  for  X  being  0  and  n-1.3  for  X  being  H  [4].  The 
number  N  of  newly  formed  X02  in  a  volume  V  is: 

N  -  V[X] [02] [M]aT-n  (2.48) 

where  [A]  denotes  the  concentration  of  the  element  A  (cm  3).  If  we 
assume  that  the  reaction  takes  place  in  an  isothermal  region  of  the 
atmosphere,  the  concentrations  can  be  represented  as: 

[X]  -  (X)/V ,  [02]  -  (02)/V ,  [M]  -  (M)/V  (2.49) 

where  the  element  symbols  in  round  brackets  denote  the  number  of 
that  element  in  volume  V.  The  equation  for  N  can  be  rewritten  as: 
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N  -  (X)(02)(M)aT_nV2 


(2.50) 


This  equation  together  with  the  thermodynamic  equation  (2.30) 
becomes : 


N 


CN 


(2.51) 


where  CN  is  a  constant  and  T  is  the  temperature.  By  incrementally 
differentiating  this  equation,  the  results  are: 


M 

N 


t-^T  "  n] 

7-1 


AT 

T 


(2.52) 


Note  that  in  the  expression  the  incrementals  AN  and  AT  are  used  in 
place  of  differentials  dN  and  dT.  This  expression  is  intentionally 
used  after  the  Krassovsky  method  in  which  the  oscillation  amplitudes 
of  the  corresponding  entities  are  referenced.  This  method  is 
discussed  in  detail  in  section  3. 

From  (2.45)  and  (2.52),  r\  can  be  represented  as: 

n  -  (2-53) 


Note  that  i;  is  a  real  number  in  this  model.  Krassovsky  gave 
the  value  of  ij  as  unity  when  X  is  an  oxygen  atom  and  3.7  when  X  is  a 
hydrogen  atom.  However,  he  reported  observed  values  of  r\  between  4 
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to  5  [Aj  .  Krassovsky  attempted  several  explanations  for  the  cause 
of  this  discrepancy,  one  of  them  being  the  need  to  consider  the 
effect  of  quenching  and  the  exclusion  of  other  chemical  reactions 
not  included  in  the  set. 

This  simple  model  has  limited  application  to  the  wave  dynamics 
of  the  mesospheric  airglow.  One  of  the  limitations  is  that  the 
measured  quantity  is  not  actually  the  volume  emission  rate  at  a 
spatial  point,  but  is  in  fact  an  integrated  value  over  a  volume 
along  the  line  of  sight  (11).  This  model  is  further  limited  by  the 
assumption  that  the  concentration  of  each  chemical  species, 
including  the  hydroxyl,  is  in  equilibrium.  This  assumption  is  valid 
only  if  the  chemical  reactions  that  determine  the  concentration 
proceed  in  a  much  shorter  time  period  compared  with  the  period  of 
the  IGW  under  investigation.  However,  several  chemical  reactions 
with  reaction  time  constants  of  5  to  20  minutes  [3]  must  be  consi¬ 
dered  when  studying  an  IGW  of  less  than  one  hour,  resulting  in  the 
introduction  of  a  phase  factor  into  rj  [11].  This  extension  of  t) 
from  an  algebraically  real  number  to  a  complex  entity  with  phase  as 
well  as  magnitude  is  discussed  in  the  following  section. 

The  work  of  Weinstock  [22],  [23]  represents  the  second  group  of 
theoretical  modeling  relating  IGWs  to  the  apparent  oscillation  in 
the  mesospheric  infrared  airglow.  Weinstock' s  model  is  a  nonlinear 
treatment  of  the  internal  gravity  wave.  This  model  is  an  attempt  to 
assess  some  apparent  nonlinear  phenomena  to  the  observations; 
however,  the  present  dissertation  does  not  include  an  implementation 
of  any  nonlinear  model. 
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The  most  obvious  difference  between  the  two  models  is  the 
dependence  of  eta  upon  the  scale  height  of  related  species, 
especially  atomic  oxygen  {11],  [23],  Hines  pointed  out  the 
difficulties  with  the  nonlinear  model  and  suggested  that  the 
apparent  nonlinearity  could  be  a  result  of  geometric  effects  [11]. 

2.  The  Complex  Eta  Model 

The  following  theoretical  calculation  of  eta  is  based  upon  that 
of  Hines  and  Tarasick  [11]  -  an  adiabatic  linearized  model  of  atmos¬ 
pheric  IGW  dynamics.  This  model  provides  a  more  detailed  con¬ 
sideration  of  the  observables  than  the  work  of  Krassovsky  (4).  This 
means  that  in  ground-based  optical  observations,  the  observed  value 
is  the  line  integral  of  the  emission  intensity  along  the  line  of 
sight  within  the  field  of  view  of  the  instrument.  As  a  result,  when 
the  line  of  sight  is  to  the  zenith,  the  horizontal  component  of  an 
incoming  wave  differs  from  the  vertical  component.  This  difference 
occurs  because  the  horizontal  (  omponent  of  the  wave  can  introduce 
additional  species  into  the  observed  column  while  the  vertical 
component  does  not.  This  is  the  same  effect  as  the  geometric  effect 
discussed  in  the  previous  section.  To  derive  tj ,  a  monochromatic, 
small-amplitude  sinusoidal  wave  of  frequency  u  is  used  as  an  input. 
Thus,  the  linearization  based  upon  the  following  formulae  can  be 
used: 


u(x,z,t)  -  u0  +  u'  exp(i(wt-kxx-kzz))  -  (2.54) 

§7  -  iuu'  (2.55) 
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- 

dx 


-ikxu' 
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(2.56) 


2u 

dz 


-ikzu' 


(2.57) 


where  u  is  a  physical  entity  such  as  velocity,  is  the  constant 
part  of  u,  u'  is  the  amplitude  of  the  sinusoidal  oscillation  of  u 
with  the  frequency  u,  the  horizontal  wavenumber  kx  and  the  vertical 
wavenumber  kz  .  The  derivation  of  the  complex  q  is  based  upon  the 
ozone  formation  reaction  (2.2): 

0  +  02  +  M  - >  03  +  M  (2.2) 

Assuming  that  the  concentration  of  the  emitting  species 
hydroxyl  is  instantaneously  proportional  to  the  concentration  of 
ozone,  the  observed  emission  intensity  is  the  integration  of  the 
volume  emission  rate  along  the  line  of  sight.  Therefore: 

00 

I  -  C  f  k*(no+n')(No+N')2  dz  (2.58) 
0 

where  k*  is  the  reaction  rate  coefficient  for  the  reaction  (2.2),  n' 
is  the  wave- induced  perturbation  of  the  concentration  [0]  from  its 
unperturbed  value  denoted  by  n0 ,  N'  is  the  wave-induced  perturbation 
of  the  concentration  of  M  from  its  unperturbed  values  No.  Here,  the 
concentration  of  02  is  assumed  to  be  proportional  to  N  at  any  time 
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with  a  ratio  C,  so  that  ( O2 ]  *■  CN.  Equation  (2.58)  can  be  rewritten 
as: 


I 


C  J  k*n0Nf  (l+jrf-)  (1+^— )2  dz 
0 


(2.59) 


The  measurable  temperature,  Tm,  will  be  taken  as  the  bright¬ 
ness-weighted  average  of  the  actual  local  temperatures  along  the 
vertical  line  of  sight.  This  can  be  represented  as: 


(2.60) 


where  T'  is  the  wave -induced  perturbation  of  T  from  its  unperturbed 

"ic 

value  T0 .  The  reaction  rate  k  is  actually  a  function  of  tempera¬ 
ture  which  can  be  linearized  as: 

k*(T)  -  k t  <J)'A 

Ao 

=  k*  (1  -  A  ■£")  (2.61) 

A0 

where  k^'  is  the  rate  constant  at  temperature  T0  and  A  is  a  constant. 
With  the  use  of  equation  (2.61),  equation  (2.59)  can  be  linearized 
and  normalized  to  10  as: 
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I-Io 
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II 

Io 


lo 


-  4-;  fc?n0HS<£-  + 


o  0 


jP  -  “4  * 

41 0  Ao 


(2.62) 


where 

00 

I0  -  /  k*n0N|  dz 

0 


(2.63) 


is  the  unperturbed  observed  intensity.  The  unperturbed  observable 
temperature  is: 


f  k*n0N|T0  dz 

To  _  0_ -  (2.64) 

f  k*n0N|  dz 
0 


Note  that  the  T0  in  the  integrand  of  the  numerator  is  a  function  of 
z..  The  fractional  perturbation  of  the  measured  temperature  is,  by 


(2.60)  and  (2.61), 


T  — X° 

m  m 


TW 

m 

TS 


S  k;n0N2(T0-TJ(El  +  -  A^)  dz 


J  ko n0NoT0  dz 
0 


(2.65) 
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where 

J  k;n0N|T'  dz 

i -  (2.66) 

f  kj  n0  N|  T0  dz 
0 

The  quotient  of  the  integrals  shown  explicitly  in  (2.65)  is  the 
terra  that  causes  the  observed  temperature  an  apparent  nonlinear 
behavior  when  the  unperturbed  temperature  T0  is  not  a  constant  with 
respect  to  altitude  z.  This  is  another  way  to  represent  the 
geometric  effect. 

The  next  step  in  the  derivation  of  rj  is  the  linearization  of 
the  continuity  equation,  here  represented  for  the  specimen  n: 


3n 

at 


-nV«v  -  y«Vn 


(2.67) 


where  v  is  the  wave-induced  velocity  vector  of  an  air  parcel,  • 
denoting  the  inner  product  of  two  vectors.  This  can  be  linearized 
with  the  assumed  input  form  of  (2.54)  into: 


(2.68) 


where  w  is  the  z  component  of  v  and  Hn  is  the  scale  height  of  the  n 
species  that  is  defined  as  follows: 
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d(ln  n0)/dz 


H 


n 


When  n  is  atomic  oxygen,  for  instance,  Hn  is  about  -3  km,  just 
under  the  OH  airglow  layer  [11].  Because  the  observables  are 
integrals  of  appropriate  physical  entities  along  the  line  of  sight, 
the  Lagrangian  method  is  more  suitable  than  the  Eulerian  [11].  In 
other  words,  each  entity  is  referred  to  with  respect  to  its  own 
stationary  coordinates  and  values.  For  example,  consider  an  air 
parcel  currently  at  a  height  z.  It  has  moved  from  its  unperturbed 
rest  level  at  height  f  -  z-h  with  the  height  increment  of  h-h(t). 
If  we  assume  that  all  species  move  together,  then  the  ratio 
n(z)/N(z)  is  the  same  as  the  unperturbed  ratio  no  (D/No  (O.  This 

gives: 

n°(0  /o  7r\\ 

n0(z)  +  n'(z)  -  (Mzj+N'^Dj^  (2'70) 

If  the  linearization  is  appropriate  in  modeling  the  travelling 
wave,  the  same  is  true  for  the  dominant  species,  namely, 

dN0(D 

N0(z)+N'(z)  -  N0(D+h"'“—  +N'(D  (2,71) 

where  any  difference  between  N' (z)  and  N’(f)  is  of  second  order  at 
most.  Since  w-iwh  by  (2.54),  (2.71)  can  be  simplified  as: 


47 


N0(z)  +  N'(z) 


N0(1+s7«v)  |j. 


(2.72) 


where  the  whole  quantity  of  the  right-hand  side  is  evaluated  at  the 
height  z-f . 

A  similar  form  for  the  kinetic  temperature  is: 

T0(z)+T'(z)  -  T0(l-4(7-l)V.v)|f  (2.73) 


where  the  adiabatic  state  and  mass  continuity  equations  were  used. 
The  detailed  linearization  of  these  equations  is  given  i.i  Appendix 
C.  Using  these  results,  the  rate  coefficient  k*  can  be  represented 
as: 


k*(z)  -  k*(z)+k'(z)  -  kfcn+h415^11  +k'(f) 

-  k}«)+[h^*r<n)[^l  |T.To(0 

-  k*[l  -  i-A(7-l)V.v]|r  (2.74) 

Then  from  (2.58)  with  (2.74), 


I0  +  I'  -  f  k*(z)[n0(z)+n'(z)](N0(z)+N'(z)]2  dz 
0 
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-  /  k*[l-iA(7-l)7.v]N3[l+iv.v]3[^-]  dz  (2.75) 

0  o  1 » 


From  the  relation: 

dz  -  (^)df  -  (l+|j?)dr  -  (2.76) 


then  (2.75)  can  be  rewritten  as  the  integration  over  f  as: 


CO 

I0  +  I'  -  f  kJn0N2[l-4(3-A(r-l))V.v-i  £*]  df  (2.77) 


which  yields: 


II 

h 


1 


J  k*n0N^[{3-A(7-l))V.y-^]  df 

CO 

J*  k*n0No  dz 
0 


(2.78) 


A  similar  linear  form  can  be  calculated  for  the  temperature 
according  to  the  following  procedure.  With  (2.74),  the  numerator  of 
(2.60)  is: 


CO 
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and  the  denominator  is  the  same  as  that  of  (2.78).  As  a  conse¬ 
quence,  equation  (2.60)  can  be  represented  in  a  normalized  form 
with  two  parts  as: 


fro 


(2.80) 


where 


s 


CO 

i  £  k*n0N|To(7-l)^«v  df 

CO 

S  kXN0T0  d2 

0 


(2.81) 


i  |  k*n0N2(T0-T«){3-A(7-l)V.v-^)  df 

CO 

J  k*n0N|T0  dz 


(2.82) 


The  first  part,  (2.81),  is  the  only  contributor  to  Tro ,  or  in  other 
words,  to  T'/T0,  when  T0  is  a  constant  for  all  z.  When  T0  is 
allowed  to  vary  with  altitude,  the  second  part,  (2.82),  works  as  an 
additional  contributor  to  Tm .  Equation  (2.81)  is  the  contributor  in 
both  cases. 

The  final  step  in  deriving  a  theoretical  co.mplex  parameter  r]  is 
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to  introduce  a  linearized  relationship  between  V*v  and  3w/3z.  Hines 
gives  such  a  relationship  [11],  [13],  namely, 

u  -  kxC2(kz-  i-(1~7^2)g]  (2.83) 

c 

w  =  w2  -  k2C2  (2.84) 

for  a  wave  of  the  form  exp[z/2H  +  i(ut-kxx-kzz) ] ,  where  g  is  the 
gravitational  acceleration,  C  is  the  speed  of  sound  and  H  is  the 
scale  height  of  the  dominant  species,  H-G~/7g-  The  horizontal 
coordinate  is  chosen  so  that  the  y  component  of  v  vanishes. 
Incorporation  of  these  relations  results  in  the  relation: 

3a  _  (k*2c2V),(-ikz+l/?H),  _ 

dz  kjj  g+ikz  u2  -w2 /2H  - 


(/i+ij/)  V»v 


(2.85) 


where  /i  and  u  are  real  and,  under  the  assumptions  made  in  this 
model,  is  independent  of  altitude.  With  the  inclusion  of  (2.85)  in 
(2.78), 


00 

i  £  k*n0N|[3-p-i^-A(7-l)]V.v  df 
°  |  k*n0N2  dz 


(2.86) 
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For  the  case  of  T0  constant  with  z,  (2.81)  and  (2.86)  together 
produce : 


H 

h 


7-1 


(2.87) 


-  ('jR+i-'jf] 


e 


(2.88) 


Equation  (2.88)  is  the  desired  complex  parameter  r)  -  r^+iqj;. 

In  the  derivation  of  (2.88),  it  was  assumed  that  the  unper¬ 
turbed  temperature  T0  is  constant  for  all  altitudes.  This  assump¬ 
tion  is  not  valid  in  the  OH  airglow  layer  region  [1].  Hines  and 
Tarasick  introduced  a  modification  of  r?  for  the  case  of  slowly 
changing  T0  with  altitude  (llj.  His  calculations  were  based  upon 
the  WKB  approximation  of  the  wave  in  a  varying  medium.  This 
approach  requires  the  following  condition: 


qH  »  1 


(2.89) 


where  q  is  equivalent  of  the  wavenumber  kz  under  the  condition  of 
spatially  variable  temperature,  and  H  is  the  scale  height.  Unfortu¬ 
nately,  this  condition  is  satisfied  by  waves  which  are  likely  to  be 
unobserved  by  ground-based  passive  optical  instruments. 

3.  Complications  of  Ela 

In  deriving  the  complex  r\  in  equation  (2.88),  it  was  assumed 
that  the  unperturbed  temperature,  T0 ,  is  constant  with  respect  to 
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altitude  in  the  region  of  interest.  In  a  better  model,  T0  should  be 
allowed  to  vary  with  altitude.  This  can  be  accomplished,  as  was 
pointed  out  by  Hines  and  Tarasick  [11],  with  the  use  of  a  WKB  appro¬ 
ximation  to  the  wave;  however,  even  with  this  improvement,  only 
those  waves  that  are  unlikely  to  be  observed  by  the  ground-based 
optical  instrument  satisfy  the  condition  (2.89)  required  by  the  WKB 
approximation  [11].  As  is  shown  by  (2.89),  those  fluctuations  which 
satisfy  (2.89)  are  the  waves  whose  vertical  wavelengths  are  much 
shorter  than  the  scale  height.  Hence,  these  waves  are  subject  to 
serious  cancellation  in  the  ground-based  passive  optical  observa¬ 
tions,  since  the  layer  thickness  is  on  the  order  of  the  scale 
height . 

The  use  of  rj  involves  another  problem.  Fig.  2-9  contains 
sample  data  taken  by  optical  equipment  described  in  Chapter  IV. 
These  airglow  measurements  clearly  show  oscillation.  Yet  from  the 
figure,  it  is  not  possible  to  confidently  ascertain  either  I0  or  I. 
The  situation  is  even  worse  for  the  rotational  temperatures  T0  or 
T'  . 

The  difficulties  described  above  should  not  be  ascribed  to  the 
definition  of  the  eta  as  given  in  equation  (2.45),  but  are  at¬ 
tributable  to  the  interpretation  of  each  element  within  the 
equation.  From  the  calculation  of  the  values  of  17  leading  to  the 
equation  (2.53),  AI  and  AT  in  (2.45)  are  incremental  fluctuations  of 
the  emission  intensity  and  the  temperature,  and  I  and  T  are  their 
absolute  values  just  before  the  fluctuations  take  place.  In  other 
words,  AI  and  AT  can  be  treated  as  differentials  dl  and  dT. 
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Instead,  in  previous  studies,  they  were  interpreted  as  the  amplitu¬ 
des  of  the  assumed  sinusoidal  oscillations  of  the  corresponding 
quantities.  By  interpreting  them  differently,  data  could  be 
processed  without  clear  definition  of  AI ,  AT,  I0  and  T0 ,  as 
discussed  in  detail  in  Chapter  VI. 

In  case  of  discrete  data,  these  derivatives  are  approximated 
by  the  increments  from  data  point  to  data  point.  A  similar  argument 
holds  for  I  and  T,  indicating  that  the  new  interpretation  is 
naturally  compatible  with  the  discrete  nature  of  the  data. 

On  the  ocher  hand,  the  departure  from  the  assumption  of  a 
sinusoidal  input  on  constant  average  values  is  undesirable.  It  was 
this  assumption  that  caused  the  cancellation  of  the  geometric 
effects  in  both  the  emission  intensity  and  the  observed  temperature 
by  the  relationship  between  V»v  and  w,  that  is,  equation  (2.85). 

One  way  to  treat  the  nonsinusoidal  input  without  constant  back¬ 
ground  is  a  direct  numerical  calculation.  The  chemical  model 
established  in  the  section  B  in  Chapter  II  can  be  used  as  a  part  of 
such  a  model.  A  numerical  model  for  the  internal  gravity  wave  is 
discussed  in  Chapter  III.  More  detailed  treatment  of  the  geometric 
effect  is  discussed  in  Chapter  IV.  Conditions  required  for  ex¬ 
perimental  data  to  be  useful  for  the  study  of  the  emission  layer  are 
discussed  in  Chapter  V. 


54 


55 


CHAPTER  III 

DYNAMIC  MODEL 


A  numerical  dynamic  model  based  upon  the  analytical  models 
developed  in  Chapter  II  is  presented  in  this  chapter.  These  methods 
include  the  use  of  MSIS  model  atmosphere  [21]  used  in  some  mesospheric 
studies  [11],  [15].  The  chemical-dynamic  model,  developed  in  Chapter 
II,  is  used  in  conjunction  with  the  MSIS  atmospheric  model  to 
calculate  infrared  emission  from  the  mesospheric  airglow.  The 
internal  gravity  wave  is  used  as  the  input  to  the  model.  Hines' 
method  introduced  in  Chapter  II  is  used  in  this  chemical -dynamic 
model.  Long-range  dynamical  characteristics  of  gravity  waves  are 
treated  in  order  to  derive  the  expected  observables.  A  model  which 
demonstrates  that  the  often  assumed  monochromatic  sinusoidal  wave  is 
not  likely  to  be  observed  in  the  mesospheric  airglow  emission  is 
examined  for  this  purpose.  This  leads  to  the  development  of  an 
optimal  signal  processing  system,  which  is  the  subject  of  Chapter  V. 

Computational  Dynamic  Model 

The  chemical -dynamic  model  introduced  in  Chapter  II  can  be  used 
with  an  atmospheric  model  that  provides  the  variable  concentrations  of 
several  species  which  form  an  emission  layer  model.  One  such  model 
can  use  the  Mass  Spectrometer  and  Incoherent  Scatter  (MSIS)  model  [21] 
which  offers  stationary  concentrations  of  N2,  03,  O,  H  and  the  average 
mass  density  p,  as  functions  of  altitude.  MSIS  also  provides  the 
temperature  for  each  altitude.  Although  MSIS  was  designed  for  therno- 
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spheric  modeling  and  may  not  be  Che  ideal  model  for  Che  mesospheric 
sCudy,  Che  following  example  shows  chac  MS1S  produces  values  close  Co 
Che  values  of  Che  US  Scandard  ACmosphere  [1].  MSIS  has  been  used  in 
several  mesospheric  sCudies  [6],  [11],  providing  one  common  reference 
for  a  comparison  of  Chese  models.  Fig.  3-1  shows  some  examples  in 
which  MSIS  produces  Che  desired  oucpuc  for  a  specific  dace. 

Using  Che  ouCpuC  from  a  model  such  as  MSIS  which  gives  che 
sCacionary  quancicies  for  Che  concenCraCions  of  chemical  species  and 
Che  Cemperacure,  Che  nexc  seep  is  Co  calculaCe  Che  concencracion  of 
Che  hydroxyl  using  Che  chemical  model  developed  in  Chapcer  II.  Fig. 
3-2  demonscraces  results  of  several  such  calculaCions .  The  alcicude 
and  Che  shape  of  Che  layer  are  consisCenC  wich  observacional  daca 
shown  in  Fig.  3-3  [10],  [24] . 

The  nexc  seep  of  Che  modeling  is  calculacing  che  dynamic  effects 
of  the  internal  gravity  wave.  As  a  first-order  approximation,  che 
wave  can  only  be  assumed  to  change  Che  concentrations  of  che  major 
species,  N2  and  O2.  As  a  second  order  approximation,  the  adiabatic 
change  of  state  due  to  IGW  is  also  included.  Once  the  change  of  state 
is  calculated,  the  chemical  model  can  be  used  to  recalculate  the 
concentration  of  the  hydroxyl  under  the  new  condition. 

In  this  study,  a  fully  detailed  computation  of  all  equations  is 
not  always  necessary  to  simulate  the  layer  emission  intensity  and 
observational  rotational  temperature.  The  reduction  method  of 
measured  airglow  data  is  given  in  Chapter  V.  It  is  necessary  to 
develop  means  of  approximation  based  upon  the  analytical  model 
developed  by  Hines,  and  introduced  in  Chapter  II  [11].  In  this  model 
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Fig.  3-3.  Measured  airglow  layer  profile  (a)  [10] 
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it  is  assumed  Chat  the  observation  is  vertical,  that  is,  the  line  of 
sight  is  toward  the  zenith.  Further  consideration  for  the  geometric 
effect  is  given  in  detail  in  Chapter  IV.  The  model  developed  in  this 
chapter  can  be  easily  modified  to  include  the  geometric  effect  for 
off-zenith  observations.  Such  an  extended  model  is  useful  for 
investigating  airglow  data  when  using  infrared  imagers  which  favor 
off-zenith  observations,  using  the  optical  brightness  gain  for 
optically  thin  airglow  layers  (5). 

Approximation  of  Fluid  Dynamics 

Among  the  three  major  steps  in  calculating  infrared  emissions, 
fluid  dynamics  is  the  most  difficult.  Calculating  fluid  dynamics 
demonstrates  the  effect  of  internal  gravity  waves  on  the  chemical- 
dynamic  atmospheric  model  by  calculating  the  temporal  changes  in  the 
concentration  of  major  species  and  also  in  the  temperature  over  small 
volumes.  This  approach  solves  three  equations  simultaneously.  These 
equations  are: 


—  -  -pVv  -  vVp 

St 


(3.1) 


3v 

p  ( —  +  (yV)v) 

a  t 


-Vp  +  pg 


(3.2) 


Vp  _  Vp 
p  ”  7  p 


(3.3) 
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where  p  is  the  mass  density,  p  is  pressure,  v  is  wind  velocity,  g  is 
the  acceleration  of  the  earth's  gravity  at  the  altitude  of  interest 
and  7  is  the  specific  heat  ratio  of  the  air  as  is  defined  in  Chapter 
II.  The  underscored  entity  is  a  vector.  These  equations  represent, 
respectively,  the  continuity  equation,  the  momentum  equation  and  the 
(adiabatic)  state  equation.  The  fluid  dynamic  equation  set,  or  the 
three  simultaneous  equations,  has  not  been  solved  in  general  form.  In 
one-dimensional  case,  some  properties  of  the  solution  have  been 
studied  (25]  and  numerical  calculation  methods  for  general  applica¬ 
tions  have  been  developed  [26],  [27].  In  two-dimensional  case,  a 
numerical  method  for  general  use  is  not  readily  available  [28]  . 
Internal  gravity  waves,  as  the  name  suggests,  naturally  involve  two 
particular  spatial  directions:  the  propagation  direction  and  the 
direction  of  the  gravitational  force.  Unlike  the  acoustic  wave,  the 
gravity  wave  treatment  employs  non- isotropic  geometry  in  which  there 
is  one  particular  direction  -  the  direction  of  the  gravity  force- 
that  renders  space  non- isotropic.  In  order  to  numerically  simulate 
internal  gravity  waves,  the  model  must  be  at  least  two-dimensional. 
The  consequence  is  that  none  of  the  existing  numerical  algorithms  for 
the  one -dimensional  case  are  applicable.  Development  of  a  two- 
dimensional  numerical  fluid- dynamics  model  is  a  subject  requiring 
intensive  study  beyond  the  scope  of  this  dissertation.  Therefore,  the 
analytical  model  developed  in  Chapter  II  is  applied  to  numerical 
calculations  in  this  chapter. 

Using  this  method,  the  original  partial-differential  equations, 
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(3.1),  (3.2)  and  (3.3),  are  converted  into  a  set  of  algebraic 
equations  with  the  assumption  that  the  input  wave  is  a  small-amplitu¬ 
de,  monochromatic,  sinusoidal  signal.  This  was  completed  in  Chapter 
II  adding  the  assumption  that  only  a  single  chemical  reaction  is 
responsible  for  the  formulation  of  the  emitting  layer.  The  results 
are  cited  as  equations  (2.80),  (2.81),  (2.82),  (2.86)  and  (2.88).  In 
this  chapter,  the  integrations  in  those  equations  are  accomplished 
numerically  using  a  more  realistic  stationary  model  generated  by  MSIS. 
The  equations  are  numerically  calculated  under  discrete  temporal  and 
spatial  steps.  The  atmosphere  is  divided  into  small  "cells"  that 
represent  local  physical  entities  such  as  temperature,  pressure  and 
species  concentration.  They  are  represented  by  (2.72),  (2.73)  and 
(2.74).  In  Chapter  II,  they  were  used  to  introduce  (2.78),  (2.81)  and 
(2.82).  Instead,  (2.72),  (2.73)  and  (2.74)  can  be  used  to  perform  the 
equivalent  numerical  integration.  This  is  the  method  used  in  this 
chapter. 

Although  it  is  not  necessary  to  carry  out  the  full  detailed 
numerical  calculations  in  order  to  study  mesospheric  airglow-emission 
oscillations,  some  consideration  which  is  related  to  the  full  scale 
model  and  necessary  for  the  adopted  simple  model  is  presented  here. 
This  involves  tha  empirical  determination  of  the  space  and  the  time 
step-size  of  the  cell  used  in  the  numerical  model. 

The  first  consideration  is  the  speed  of  sound  in  the  mesospheric 
region,  as  shown  in  fig.  3 -1(f).  The  speed  of  sound  is  the  speed  at 
which  the  displacement  of  pressure  from  its  stationary  value  propaga¬ 
tes.  This  is  the  fastest  speed  at  which  any  disturbance  propagates. 
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Therefore,  Che  ratio  of  the  spatial  step  to  the  time  step  should  be 
kept  on  the  order  of  the  value  of  the  sound  speed. 

The  next  consideration  is  the  wind  speed  in  the  region.  Many 
observational  data  record  a  wind  speed  of  as  much  as  50  to  100  m/s 
around  the  emitting  layer  altitude  [5],  [13].  Therefore,  the  spatial 
step  can  be  several  hundred  meters  to  a  few  kilometers. 

The  third  consideration  is  the  typical  oscillation  period  of  an 
internal  gravity  wave.  Theoretically,  as  explained  in  Chapter  II,  the 
period  of  a  linearized  atmospheric  IGW  is  longer  than  the  Brunt - 
Vaisala  period,  which  is  about  5  minutes  in  the  mesospheric  region,  as 
shown  in  Fig.  3-l(g).  Consequently,  the  time-step  could  be  as  long 
as  several  minutes.  With  these  considerations,  the  spatial  step  of 
one  kilometer  and  the  time  step  of  5  to  10  seconds  could  have  been 
selected.  However,  as  Fig.  3-1  shows,  the  values  of  most  of  the 
stationary  background  quantities  can  change  significantly  in  one 
kilometer.  Also,  as  Fig.  3-2  and  3-3  indicate,  the  infrared  emission 
layer  is  about  7  km  thick;  therefore,  a  spatial  step  of  much  less  than 
one  kilometer,  such  as  100  meters,  is  a  more  desirable  choice.  This 
causes  the  corresponding  time  step  to  be  about  one  second.  Note  that 
this  process  of  cell-size  determination  is  also  valid  for  full-scale 
numerical  simulation  of  the  original  differential  equations,  (3.1), 
(3.2)  and  (3.3). 

In  the  case  of  small-amplitude  IGW's  observed  in  the  mesosphere, 
the  wave  can  be  considered  quasi-monochromatic.  The  justification  is 
that  such  waves  are  believed  to  be  generated  at  low  altitudes  and  are 
then  propagated  through  the  stratosphere  and  lower  mesosphere  [13]. 
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In  these  regions  of  IGW  propagation,  there  exist  temperature  inver¬ 
sions.  There  is  also  a  dispersion  relation  between  the  wave  period 
and  the  propagation  direction  [18J.  These  two  conditions  result  in  a 
"window"  for  the  wave  frequencies  allowed  in  the  mesopause.  As  a 
result,  an  IGW  generated  at  low  altitude  is  expected  to  produce  close 
to  a  monochromatic  wave  in  the  upper  mesosphere.  This  justifies  the 
assumption  that  locally  monochromatic  waves  exist  in  the  mesosphere. 
Note,  however,  that  this  assumption  is  valid  only  when  the  amplitude 
of  the  wave  is  small  enough  so  that  the  linearization  of  the  fluid- 
dynamic  equations  are  valid.  When  the  wave  has  a  large  amplitude,  the 
lower  atmosphere  may  no  longer  act  as  a  band-pass  filter,  possibly 
allowing  non-monochromatic ,  or  even  non-sinusoidal ,  waves  to  reach  the 
mesosphere.  The  linear  model  is  inadequate  to  describe  such  waves  of 
great  amplitude.  Such  cases  need  a  nonlinear  treatment  which  is 
beyond  the  scope  of  this  dissertation.  As  pointed  out  in  Chapter  VII, 
solitary  waves,  or  "solitons",  and  a  mathematical  chaos  model  are 
among  possibilities  to  be  studied  in  future.  IGW  dynamics  leading  to 
quasi -monochromasy  in  the  mesosphere  are  discussed  in  subsequent 
sections . 

Approximation  of  Chemical  Model 

The  chemical  model  equations  (2. 4) -(2. 6)  described  in  Chapter  II 
may  be  used  to  perform  a  numerical  calculation  using  sufficiently 
small  time -step  At.  Obviously,  the  At  selected  must  be  much  shorter 
than  the  time  constant  of  the  fastest  chemical  reaction  involved  in 
this  calculation.  The  fastest  chemical  reaction,  represented  by 
equation  (2.5),  has  a  time  constant  of  tens  of  milliseconds.  Conse- 
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quently,  At  was  selected  to  be  less  than  or  equal  to  one  millisecond. 
With  a  time  step  of  one  millisecond,  it  takes  one  million  calculation 
steps  to  simulate  ten  minutes  in  real  time. 

Another  method  is  available  to  calculate  the  chemical  reactions 
(2.4)-(2.6),  based  upon  the  following  assumptions: 

1]  Diffusion  of  all  species  is  insignificant. 

2]  The  concentrations  of  the  major  species  remain  almost  constant 
within  the  IGW  period. 

3]  The  chemical  reactions  which  determine  the  concentrations  of  the 
active  minor  speci.es  take  place  in  a  time  scale  much  shorter 
than  the  time  scale  of  interest.  Therefore,  the  concentrations 
can  be  regarded  as  in  equilibrium. 

In  other  words,  the  effects  of  the  transient  behavior  of  the 
active  minor  species  are  assumed  to  be  negligible  within  the  IGW 
period.  Based  on  these  assumptions,  the  concentrations  of  the 
chemically  active  minor  species  of  interest  are  represented  as 
functions  of  the  concentrations  of  the  major  species.  This  form  puts 
forth  the  desired  concentrations  in  terms  of  stable  and/or  measurable 
values.  Therefore,  a  "hybrid"  approach,  the  combination  of  the  direct 
numerical  simulation  and  the  equilibrium  value  model,  is  feasible. 
When  combined  with  the  dynamic  numerical  model  introduced  previously, 
the  equilibrium-value  approach  significantly  reduces  the  amount  of  the 
total  calculation.  Furthermore,  if  the  concentration  of  the  hydroxyl 
alone  is  assumed  to  be  in  equilibrium,  the  next  fastest  reaction  is 
that  of  ozone,  (2.3),  with  a  time  constant  of  tens  of  minutes.  This 
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means  that  even  if  the  '•  nsient  behavior  of  the  ozone  concentration 
is  to  be  calculated,  the  time  step  may  be  set  as  long  as  tens  of 
seconds.  As  a  consequence,  the  total  number  of  calculations  is 
reduced  by  a  factor  of  ten  thousand,  that  is,  from  one  millisecond  to 
ten  seconds  for  At.  The  At  was  selected  at  about  one  second  in  the 
dynamical  model,  a  value  which  is  also  appropriate  for  the  chemical 
model  used  in  the  above  manner. 

In  most  cases,  the  calculation  of  equilibrium  values  of  minor 
species,  assuming  that  long-lifetime  species  concentration  is  quasi¬ 
constant,  is  sufficiently  accurate  compared  to  the  direct  numerical 
calculation.  However,  it  fails  to  produce  accurate  results  in  cases 
where  the  system  is  not  in  equilibrium,  such  as  at  twilight  [3],  In 
other  words,  this  hybrid  approach  should  be  used  when  the  time  period 
of  interest  is  about  30  minutes  or  longer.  When  the  time  period  of 
interest  becomes  shorter  than  about  30  minutes,  equation  (2.3)  should 
be  numerically  calculated  for  each  At. 

Geometric  Effect 

The  geometric  effect  plays  a  very  important  role  in  modeling  of 
the  observables  of  the  mesospheric  airglow.  To  include  this  effect  in 
the  numerical  models  introduced  thus  far,  the  most  straight-forward 
method  is  the  numerical  integration  of  appropriate  quantities  along 
the  line  of  sigh*-.  Since  the  chosen  spatial  step  of  100  meters  is 
small  enough  to  assume  that  most  relevant  physical  quantities  are 
constant  within  one  cell,  the  approach  can  be  -further  simplified. 
Instead  of  numerical  integrations,  appropriate  quantities  are  summed, 
aus-  •  ".  has  been  found  (Chapter  II)  that  the  geometric  effect  has 
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significant  influences  on  many  observables,  a  more  detailed  analysis 
of  the  effect  is  needed.  Because  off-zenith  angle  observations  are 
often  exercised,  the  geometric  effect  of  oblique  observations  should 
also  be  examined.  These  considerations  are  detailed  in  the  next 
chapter . 

Model  Calculation  Examples 

In  this  section,  several  results  of  calculations  utilizing  the 
hybrid  model  introduced  in  the  previous  section  are  presented.  The 
background  physical  quantities  obtained  by  the  MSIS  model  are  used  and 
the  linearized  equations  (2.80),  (2.81),  (2.82)  and  (2.86)  are 
numerically  calculated.  Outputs  consist  of  integrated  infrared 
emissions  and  intensity-weighted  temperature  along  the  vertical  line 
of  sight.  The  geometric  effect  for  the  vertical  observation, 
therefore,  is  already  taken  into  account.  Note  that  the  discussion  of 
r)  which  led  to  the  equations  (2.87)  and  (2.88)  requires  no  numerical 
calculation  with  respect  to  time,  since  the  integration  involving  the 
exp(iwt)  term  is  cancelled  and  does  not  appear  in  the  latter  equa¬ 
tions.  Calculations  reported  here  demonstrate  typical  signal  forms 
that  are  expected  to  be  observed  according  to  the  numerical  model. 

The  integration  is  implemented  by  adding  the  real  part  of  the 
complex  integrands  in  (2.80),  (2.81)  and  (2.86),  since  the  real  part 
represents  the  physical  entity.  The  phase  term  exp  i(wt-kxx-kzz) 
changes  as  time  progresses,  and,  for  the  vertical  line  of  sight,  x  is 
constant  for  all  z  so  that,  without  a  loss  of  generality,  x  can  be  set 
to  0.  This  term  is  now  simply; 
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V»V  -  (cos(wt-kzz)  +  i  sin(wt-kzz))  (3.4) 

This  equation  plays  a  major  role  when  off-zenith  observations 
are  simulated  using  the  same  numerical  model.  This  is  explained  in 
detail  in  Chapter  IV. 

1  _g.  3-4  presents  several  results.  These  examples  demonstrate 
the  geometric  effect  on  the  intensity  fluctuation.  The  thickness  of 
the  hydroxyl  layer  is  determined  by  the  scale  height  of  related 
chemical  species,  which  is  a  function  of  temperature.  The  ratio  of 
the  layer's  thickness  to  the  vertical  wavelength  of  an  IGW  dictates 
the  modulation  occurring  in  the  integrated  emission  intensity.  The 
modulation  in  temperature,  however,  is  not  affected  in  the  same 
manner.  The  temperature  fluctuation  is  influenced  by  the  value  of  the 
stationary  temperature,  not  by  the  ratio  of  the  layer  thickness  to  the 
IGW  wavelength.  This  example  demonstrates  the  geometric  effect 
introduced  in  Chapter  II  under  more  realistic  conditions.  However, 
the  result  is  consistent  with  the  simplified  model  of  Chapter  II. 

Long  Range  Dynamics  of  IGW 

Because  the  origins  of  IGW's  are  believed  to  be  distant  from  the 
mesospheric  region  where  the  observational  data  are  taken,  it  is 
important  to  examine  the  characteristics  of  the  excitation  and  the 
long  range  propagation  of  IGWs.  f  ..ch  an  analysis  will  assist  in 
predicting  signal  characteristics  for  ground-based  optical  observa¬ 
tions.  The  models  developed  so  far  use  monochromatic  sinusoidal  waves 
as  input.  It  will  be  demonstrated  in  this  section  that  waves  found  in 
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Fig.  3-4.  Numerical  model  calculation 
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Che  mesospheric  region  are  not  likely  to  be  monochromatic,  unless  the 
source  of  the  IGW  is  a  monochromatic  sinusoidal  oscillator.  The 
signal  characteristics  are  determined  by  the  dynamics  of  the  IGW 
generation  and  propagation,  as  well  as  by  the  geometrical  configura¬ 
tion  of  the  observation. 

In  this  section,  the  effect  of  the  dynamic  propagation  of  an  IGW 
is  discussed  by  showing  that  the  expected  signal  is  chirp  sinusoidal 
as  a  result  of  these  mechanisms.  Details  of  the  geometric  effect  are 
left  to  Chapter  IV. 

The  propagation  of  linearized  IGVs  summarized  by  Hines  (18}  is 
presented  in  Appendix  A.  Only  the  results  are  used  here.  In  this 
model,  the  source  of  IGW  is  assumed  to  be  in  lower  atmosphere  -  the 
upper  troposphere  or  the  lower  stratosphere.  In  many  cases,  this 
source  appears  in  the  form  of  an  atmospheric  discontinuity,  such  as  a 
front  system,  a  (thunder)  storm  or  strong  winds  colliding  with  tall 
mountains  (5).  Once  an  atmospheric  disturbance  is  generated,  its 
energy  propagates  through  the  atmosphere  at  the  packet  velocity 
(Appendix  A)  .  In  the  course  of  propagation,  the  dispersion  relation 
represented  by  equation  (A. 4)  in  Appendix  A  causes  the  wave  packet  to 
spread  spatially.  The  separation  occurs  in  the  propagation  packet 
speed  according  to  the  propagation  direction  and  the  oscillation 
frequency.  Since  the  hydroxyl  layer  at  the  stationary  state  can  be 
modeled  to  be  flat  at  a  constant  altitude,  the  observable  of  ground- 
based  optical  measurement  is  the  intersection  of  such  a  horizontal 
layer  with  the  propagating  wave.  Note,  however,  that  such  an 
interaction  can  be  observed  as  the  local  phase  velocity.  The  packet 
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velocity  can  not  be  measured  using  a  ground-based  optical  instrument 
with  a  small  (1  degree  or  so)  field  of  view.  As  shown  in  Appendix  A, 
in  cases  in  which  IGWs  are  generated  at  low  altitudes  propagating  with 
an  upward  packet  velocity,  the  phase  velocity  is  downward  and  the 
phase  speed  is  greater  than  the  packet  speed.  Consequently,  the  cross 
section  of  the  airglow  layer  interacting  with  the  IGW  should  appear  as 
shown  in  Fig.  3-5.  The  slanted  phase-plane  interacting  with  an 
airglow  layer,  .  whose  thickness  is  comparable  with  the  vertical 
wavelength  Az,  causes  the  geometric  effect. 

Since  the  origins  of  many  IGWs  are  forms  of  discontinuities,  the 
generated  signal  is  likely  to  contain  broad-band  frequency  components!' 
As  depicted  in  Fig.  3-6,  each  frequency  component  has  a  unique 
"resonant  direction,"  i.e.  a  direction  towards  which  energy  propagates 
faster  than  in  other  directions.  The  equations  upon  which  the  figure 
is  based  are  introduced  in  the  following  section. 

For  each  direction,  a  particular  wave  frequency  component 

propagates  faster  than  other  frequency  components,  as  depicted  in 
Fig.  3-7  [18].  These  properties  indicate  that  after  travelling 

several  hundred  kilometers,  a  wave  packet  which  has  close  frequency 
components  is  likely  to  develop.  The  center  frequency  of  each  packet 
is  a  function  of  the  direction  of  propagation. 

Such  wave  frequency  components  are  subject  to  the  filtering 

effect  of  the  temperature  inversion,  as  described  in  Chapter  II. 

Tims ,  only  certain  frequency  components  propagate  through  a  tempera 
ture  Inversion  and  ultimately  reach  the  airglow  layer.  Fluctuations 
produced  by  such  wavepackets  in  the  airglow  layer  are  sinusoidal -like 
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Fig,  3-6.  Packet  velocity  vs.  elevation  angle 
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oscillations  with  slowly  changing  frequency,  often  referred  to  as 
chirp  signals  (Fig.  3-8).  Note  that  the  quasi -raonochroina tic  wave 
model  can  be  still  used  cautiously  in  many  cases,  because  the 
frequency  usually  changes  only  a  few  percent  during  one  cycle.  Also, 
experimental  data  are  obtained  only  for  a  few  cycles  -  in  other  words, 
the  amount  of  observational  data  is  small  enough  to  be  treated  as 
monochromatic.  Such  slow  frequency  shifts  are  reported  from  observa¬ 
tions  [5] . 


Model  Kesul ts 

The  packet  velocity  P-(PX,PZ)  at  which  a  wave  packet  of 
frequency  w  propagates  is: 
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k2  w3 


When  this  wave  packet  is  propagating  with  the  elevation  angle  a, 
the  following  equation  is  also  satisfied: 

kz  -J  P  i 

tan  a  ~  —  -  ~  —  (A.  16) 

kx  u2  -  u»|  Px 

When  a  wave  packet  has  a  constant  propagation  angle  a,  the 
equation  (A.  14)  can  be  used  to  calculate  Px,  and  then  use  of  the 
relation  P2-Px  tana  yields  Pz .  In  order  to  calculate  Px  for  a  given 
u,  kx  needs  to  be  known.  This  can  be  accomplished  by  solving 
equations  (A. 4)  and  (A. 16),  which  yields: 

2  (u2-u2)2 

(w2-  KUg )  —  +  Wgk2  -  w2k2(  1  +  -  tan2  a)  -  0  (3.5) 

C2  J' 

where  /c-w3/wg=l . 225  (18).  For  a  given  u,  kx  is  calculated  numerically 
using  (3.5),  and  with  the  w  and  resultant  kx,  (A.  14)  can  be  calculated 
numerically.  Solving  (3.5)  for  kx  yields: 

w*  ('x>2  -  <w2) 

k2  - - - -  (3.6) 

C2  (w2-Wg)(w2  +  (u2-w2 ) tan2  a) 
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For  an  IGW,  w  <  ug  .  Therefore,  for  kx  to  be  algebraically  real, 
w2  +  (w2_w2 )  tan2a  >  0.  This  can  be  rewritten,  in  terms  of  the  wave 
period  T-2 n/u  and  Tg-27r/ug  as: 

1 

T2  <  T2  (  1  +  -  )  (3.7) 

tan2  a 

This  constraint  indicates  that  for  a  given  wave  propagation  angle 
there  is  a  maximum  possible  oscillation  period,  or  minimum  frequency, 
for  physically  allowed  IGW.  Fig.  3-9  plots  the  maximum  period  given 


Tmax 
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Fig.  3-9.  Limiting  period  and  elevation  angle 
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by  the  right  hand  side  of  equation  (3.7),  Tmax,  versus  a;  the  same 

function  is  plotted  in  several  different  scales.  For  each  a,  any 

wave  component  with  a  period  smaller  than  T  can  propagate.  Any 

max 

wave  component  with  a  period  larger  than  Tmax  can  not  propagate  at  the 
angle  a.  In  other  words,  all  frequency  components  have  upper  limits 
of  elevation  angle  a,  represented  in  the  right  hand  side  of  equation 
(3.7). 

Fig.  3-10  shows  several  calculated  values  of  Px.  Each  frame  is 
a  plot  of  wave  packet  speed  Px  as  a  function  of  wave  frequency, 
represented  in  period  T-2tt/w.  The  ascent  angle  a  varies  from  frame  to 
frame . 

The  time  Tx  for  a  wavepacket  with  speed  Px  to  reach  at  distance 
Lx  is  Lx/Px.  The  observation  of  an  IGW  takes  place  at  a  layer  of 
altitude  h.  Using  the  right  hand  side  of  equation  (3.7),  the  distance 
between  the  wave  source  and  the  observation  location  can  be  calculated 
by  Lx-h/tan  a.  Table  VIII  is  an  example  of  a  and  calculated  Lx  . 

Fig.  3-11  is  a  plot  of  the  frequency  variation  resulting  from 
the  interaction  between  an  IGW  and  a  horizontal  layer  at  an  altitude 
of  85  km.  One  important  feature  is  that  at  each  moment,  there  are  two 
w  values.  This  indicates  that,  except  for  the  frontal  frequency  (i.e. 
the  frequency  that  reaches  to  the  observation  location  first)  a  pair 
of  frequencies  may  be  observed. 

A-lso,  the  frequency  change  rate  is  a  function  of  frequency 
(Fig.  3-11).  For  each  a,  there  are  two  cut-off  frequencies:  the 
Brunt-Vaisala  frequency  at  the  short  period  limit  and  the  limiting 
frequency  corresponding  to  the  limiting  period  Tmax  represented  by  the 
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Fig.  3*10.  Packet  velocity  (continued) 


tn  n  CM  H 


period 


0 


60 


lime  (min) 


m 


Fig.  3-11.  Frequency  chirp  (shown  with  period)  (continued) 
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Fig.  3-11.  Frequency  chirp  (shown  with  period)  (continued) 

right  hand  side  of  equation  (3.7),  as  plotted  in  Fig.  3-9.  The 
frequency  change  rate  is  maximum  at  the  frontal  frequency  and  becomes 
moderate  at  its  neighboring  frequencies.  As  time  passes,  the  frequen¬ 
cies  become  almost  constant,  corresponding  to  low  and  high  limiting 
frequencies . 

These  results  indicate  that  an  IGW  which  was  generated  in  a  distant 

/ 

low  atmosphere  and  has  propagated  to  the  mesosphere  is  likely  to  chirp. 
It  is  therefore  necessary  to  construct  a  system  which  will  process  such 
signals.  Such  a  system  is  developed  in  Chapter  V,  and  evaluated  in 
Chapter  VI . 
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CHAPTER  IV 

GEOMETRIC  EFFECT 


In  this  chapter,  the  geometric  effect  is  studied  in  detail.  This 
effect,  introduced  in  Chapter  II  in  a  very  simplified  manner,  has  been 
shown  to  have  a  significant  effect  on  observables;  the  airglow 
emission  intensity  and  the  rotational  temperature.  Also  included  in 
this  chapter  is  the  geometric  effect  of  the  off-zenith  observation. 

Geometric  Effect  Model 

Among  the  several  ways  to  calculate  the  observables  of  interest 
using  a  model  with  geometric  effect  taken  into  account,  one  straight¬ 
forward  method  involves  the  numerical  integration  of  appropriate 
quantities  along  the  line  of  sight.  The  two  observables  of  interest 
are  the  infrared  airglow  emission  intensity,  or  the  radiance,  and  the 
rotational  temperature  (Chapter  V  and  Appendix  B) .  The  rotational 
temperature  can  be  assumed  to  be  identical  to  the  kinetic  temperature 
of  the  major  species  of  the  emission  layer  region,  if  the  time  L.cale 
of  interest  is  at  least  a  few  seconds.  When  integrated  through 
altitude,  however,  the  observed  kinetic  temperature  is  not  always 
identical  to  the  observed  rotational  temperature.  This  occurs  because 
the  rotational  temperature  is  calculated  from  the  observed  distribu¬ 
tion  of  intensities  of  the  vibrational-rotational  spectral  lines,  and 
each  spectral  line  is  integrated  through  the  line  of  sight  indepen¬ 
dently.  The  rotational- line  intensity  is  a  function  of  the  total 
emission  intensity  and  kinetic  temperature.  Since  the  airglow  layer 
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density  and  kinetic  temperature  change  independently  in  real  atmo¬ 
sphere,  distribution  of  the  integrated  line  intensity  does  not 
correspond  to  any  total  intensity  and  kinetic  temperature  combination. 
For  comparison  purpose,  both  the  observed  kinetic  temperature  and  the 
rotational  temperature  arc  calculated  according  to  the  geometric 
effect.  The  discrepancy  between  the  two  temperature  is  then  demon¬ 
strated. 

The  formulae  for  the  integration  through  altitude  are: 


J  I(z,t)  dz 
0 


(4.1) 


9  i-1 

l(z, t)  -  S  { [0H*(z) ]  S  A(i,j)l 
1*1  j*1 


(4.2) 


Ot)  -  S  T(z,t)I(z,t)  dz 


obs 


(4.3) 


where  "kin"  stands  for  kinetic,  "obs"  stands  for  observable,  or 
observed.  These  formulae  represent  the  same  principle  used  in  the 
introduction  of  the  geometric  effect  in  Chapter  II. 

The  method  of  extending  the  calculations  for  the  above  formulae  to 
off-zenith  observation  cases  is  based  upon  the  Hines  model  introduced 
by  equations  (2.72)  and  (2.73).  The  basic  principle  was  to  convert 
the  variable  from  the  actual  altitude  z  to  its  stationary  altitude  f . 
This  method  assumes  monochromatic  sinusoidal  gravity  wave  input.  The 
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first  step  involves  calculating  amplitudes  of  the  fluctuation  in  the 
concentration  of  the  major  species  and  the  kinetic  temperature  for 
each  f  using  the  Hines  model.  Using  those  amplitude  and  the  desired 
vertical  wavenumber,  the  temperature  and  the  concentration  for  each  f 
at  time  t  are  calculated.  This  can  be  done  by  multiplying  the  term: 

exp(i(ut-k„x~kzz))  (4.4) 

to  the  integrand  for  each  f  derived  by  the  model,  then  adding  the  real 
parts  of  the  resultant  complex  values  according  to  equations  that  are 
converted  from  (4.1),  (4,2)  and  (4.3),  changing  z  to  f  to  obtain  the 
desired  quantities.  The  equations  are  (2.80),  (2.81),  (2.82)  and 
(2.86).  The  theoretical  q  can  be  calculated  by  (2.88).  Later 
sections  demonstrate  that  this  method  allows  a  simple  extension  of 
the  Hines  model  to  a  case  of  off-zenith  observation. 

Geometric  Effect  on  Rotational  Temperature 

There  are  several  ways  to  calculate  the  rotational  temperature 
from  the  spectral  line  data.  One  of  such  methods  is  to  fit  a 
temperature  to  two  of  such  line  data  according  to  the  following 
formulae  (10] : 


C(t) 


(he)  (Ejl  -  Ej2) 

ln(IiLAil) 

Ii2  Ail 


(4.5) 
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where  Iji  is  the  rotational  line  intensity  of  the  emissive  transition 
from  (v'ji)  to  (v"j"),  Eji  is  the  energy  level  of  the  v'ji  state 
represented  in  the  unit  of  cm'^-.  is  its  Einstein  coefficient  for 
the  emissive  transition  and  jl  and  j2  are  two  different  rotational 
levels  with  the  same  upper  vibrational  level.  These  formulae  are 
based  upon  the  theoretical  spectral  form  of  {9],  {10]: 


I(NV,  .T.i/.j,)  - 


NV  ^ii 1  SH‘ 


2 ( 2 j +1 ) exp ( -E  t /kT ) 
j=1/2  1 


exp(-Ej/*T) 


(4.8) 


where  I  is  the  volume  emission  rate  (Rayleigh)  of  the  emission  line  of 
the  wavenumber  »/..(  (cm  ),  i  and  i'  denote  the  upper  and  lower 
rotational  state,  at  temperature  T  (K)  with  the  number  density  of  the 
upper  vibrational  state  hydroxyl  Nv'  (cm*^).  S^,  is  the  line 
strength  and  the  denominator  of  (4.7)  is  the  partition  function  for 
the  upper  vibrational  state.  Theoretical  values  with  the  corrections 
using  the  observed  data  should  be  used  for  parameters  v  ,  S.{1  and 
E j .  Such  data  are  available  from  several  sources  [9],  [29],  [30]. 
In  this  study,  the  data  developed  at  Utah  State  University  based  upon 
Heis  model  is  used  [29]. 
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Another  way  to  calculate  the  rotational  temperature  is  to  fit  a 
temperature  to  more  than  three  pairs  of  E^(  and  ln(I  A  )  where  j' 
represent  different  rotational  levels  of  the  same  upper  vibrational 
level  v' .  The  fitting  of  the  temperature  can  be  done  with  the  least- 
square  technique  [31],  This  section  demonstrates  that  a  large 
fitting  error  could  result  even  from  data  with  a  high  signal  to  noise 
ratio,  when  the  geometric  effect  is  taken  into  account. 

The  least-square  fitting  of  form  y-ax+b  to  the  data  (x^.y^)  with 
the  standard  deviation  error  for  (y^)  being  (eri )  (i-l,N)  has  a  form 
[32]: 


SxxSy-SxSxy 

A 


(4.9) 


b  _  SSxy-SxSy 
A 

with  the  standard  deviation  error  estimate 


Sxx 


(4.10) 


(4.11) 

(4.12) 


where 


A  -  SSxx-(Sx)2 


(4.13) 
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In  case  of  the  rotational  temperature  calculation,  the  data  pair  (x,y) 
should  be  x-hcEj/k  and  y-ln(I^/i/^.'S..*) .  The  linear  fitting  can  Chen 
provide  appropriate  temperature-error  estimates. 

Both  the  two-lines  method  and  the  least-square  fit  method  assume 
that  the  rotational  emission  lines  are  results  of  infrared  emitting 

hydroxyl  volume  of  a  single  temperature.  This  is  not  the  case  in  a 

\ 

ground-based  passive  observation,  when  there  is  a  temperature 
difference  within  the  emitting  layer.  This  is  apparent  from  (4.6)  and 
(4.8)  combined  into  the  following  formula: 


Iji  "  /  I(Nvl(z),T(z),.,fJ,)  dz 


®  N  (z)  v  S 

J  - - - - - — -  exp(-Ej/«T(z))  dz 

0  2(2j+l)exp(-E:AT(z))  J 

M/2  J 


(4.15) 


For  comparison,  both  the  brightness -weighted  temperature  and  the 
rotational  temperature  calculated  from  the  integrated  rotational  lines 
are  presented  in  the  next  section. 
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The  significance  of  the  temperature  change  in  the  emission  layer 
and  the  emission  layer  profile,  represented  by  T(z)  and  Nv. (z)  in  the 
numerator  of  (4.15),  is  demonstrated  using  the  computer  model.  This 
model  provides  the  integrated  line  intensities  for  each  rotational 
line  according  to  (4.15)  and  the  rotational  temperature  estimate 
according  to  the  least  square  fit.  The  model  also  provides  the 
intensity  weighted  kinetic  temperature  according  to  (4.3).  For  the 
emission-layer  profile,  the  results  of  the  Gaussian  distributions 
with  different  variance  above  and  below  the  peak  and  the  distribu¬ 
tions  derived  from  MSIS  are  compared.  The  kinetic  temperature 
profile,  parabolic  profile  with  different  coefficients  and  profile 
derived  by  MSIS  are  used  and  compared.  The  resultant  geometric 
effects  are  shown  quantitatively. 

In  this  model,  the  volume  emission  rate  profile  is  defined  as: 

(OH(z)  ]  -  exp(—  £3--Z^2-)  (4.16) 

a2 


for  the  simple  Gaussian  profile  with  the  peak  altitude  zO  and 
thickness  parameter  a,  or 


[0H(z)J  -  exp(-  -k-aQ>2.)  (zszO) 

(4.17) 

-  exp(-  (z>z0) 
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for  the  Gaussian  with  different  thickness  parameters  cL  and  a y ,  below 
and  above  the  peak  altitude  zO,  respectively.  Examples  shown  in  Fig. 
4-1  and  4-2  provide  close  approximations  to  the  MSIS  model  introduced 
in  Chapter  III,  shown  here  for  comparison  purposes. 

The  temperature  profile  is  defined  as: 

T(z)  -  T0  +  AT(z-zt)2  (4.18) 


for  a  simple  parabolic  profile  with  the  minimum  temperature  Tg  at 
altitude  zt  and  a  parameter  AT.  To  make  this  model  closer  to  MSIS 
model,  the  following  model  is  used: 


T(z)  -  T0  +  ATL(z-zt)2 


<sizt> 


-  T0  +  AT0  (z-zt  )2 


(zazt) 


(4.19) 


which  is  a  parabolic  profile  with  the  minimum  temperature  T0  at 
altitude  zt  and  different  parabolic  parameters  ATL  and  AT0  for  below 
and  above  zt  ,  respectively.  Examples  shown  in  Fig.  4-3  and  4-4 
provide  close  approximations  to  the  MSIS  model  presented  in  Chapter 
III  and  included  in  the  figure. 

In  these  models,  each  rotational  line  emission  intensity  for  each 
altitude  is  calculated  using  equation  (4.8).  The  resultant  inten¬ 
sities  are  added  according  to: 


I  (j,k)  -  Az  S  I(Ni,Ti,^iU) 

obs  i=0  J 


(4.20) 
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where  j  and  k  denote  rotational  levels,  N^  and  denote  the  con¬ 
centration  of  hydroxyl  at  the  upper  vibrational  level  v*  and  tempera¬ 
ture  at  altitude  z^,  and  z^-z^+iAz.  In  the  following  examples,  z^  is  . 
set  to  70  km,  Az  is  0.1  km  and  i  runs  from  0  to  301,  resulting  z^  from 
70  to  100  km.  is,  therefore,  the  local  volume  emission 

rate  of  wavenumber  at  altitude  Z£. 

Results  of  this  model  are  provided  in  Fig.  4-5  and  Table  IX. 
Fig.  4-5  contains  plots  of  ( j ) ,  the  energy  of  the  initial  rota¬ 
tional  level  j,  and  natural  logarithm  of  I^(j ^  where 
denotes  wavenumber  and  denotes  the  line  strength  of  the  rotational 
emission  line  Table  X  lists  parameters  of  several  emission 

lines  ordered  by  their  energy.  Table  XI  contains  the  estimated 

intensity  and  temperature  with  their  error  estimates,  corresponding  to 
Fig.  4-5.  Table  XI  also  includes  the  integrated  intensity  weighted 
kinetic  temperature.  Note  that  the  rotational  and  kinetic  tempera¬ 
tures  are  close  but  not  identical  in  most  cases.  Fig.  4-6  provides 
rearrangement  of  the  emission  intensity  and  temperature.  Instead  of 
functions  of  altitude,  the  intensity  at  each  temperature  is  accumula¬ 
ted  for  all  altitude  and  plotted  the  apparent  intensity  weight  of  each 
temperature.  If  there  is  a  passive  optical  instrument  which  can 
measure  kinetic  temperature,  such  instrument  will  measure  the 


intensity  distribution  as  a  function  of  temperature,  similar  to  this 
figure.  It  is  clear  that  the  distribution  is  far  from  gaussian  in 
many  cases.  Those  cases  with  small  temperature  error  estimates 
somehow  retain  symmetry  and/or  small  variances',  while  those  with 
large  error  estimate  are  asymmetric  and/or  have  large  variances.  This 
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TABLE  IX 

Rotational  Temperature  Fit  Statistics 


j 

EiU) 

ln(IAjkSjk) 

CASE  1:  1 

0.0 

1.38540 

2 

107.477 

0.68924 

3 

259.126 

-0.29077 

4 

455.375 

-1.55459 

5 

696.371 

-3.09882 

CASE  2:  1 

0.0 

1.54621 

2 

107.477 

0.85288 

3 

259.126 

-0.12264 

4 

455.375 

-1.37978 

5 

696.371 

-2.91419 

CASE  3:  1 

0.0 

1.53226 

2 

107.477 

0.84781 

3 

259.126 

-0.11533 

4 

455.375 

-1.35669 

5 

696.371 

-2.87224 

TABLE  X 

Rotational  Line  Parameters 


111 


Energy  (cm'1) 

E0  ->  10172.28 

10172.32 
10246.99 
10247.15 

10300.34 
10300.48 
10352.27 
10352.63 
10354.09 

10354.32 
10443.17 

.  10443.42 

10488.46 

10489.11 
10566.93 

10567.11 

10655.80 
10656.82 
10724.70 
10724.73 

10854.37 

10855.84 
10915.67 

10915.88 

11084.14 

11086.14 

11139.37 

11139.89 

11344.96 
11347.57 
11395.30 
11396.22 
11636.61 

11639.89 

11682.96 

11684.35 

11958.78 

11962.81 

12001.85 

12003.78 
12351.44 
12353.98 


E-E0  (ea* 1 )  j 

0.0000000E+00  2 

4.2968750E-02  2 

74.71387  3 

74.87305  3 

128.0635  1 

123.1592  1 

179.9873  * 

180.3477  * 

181.8135  2 

182.0420  2 

270.8396  3 

271.1337  3 

316.1757  5 

316.8251  5 

394.6475  4 

394.8201  * 

483.5186  6 

4S4.54C0  6 

552.4229  5 

552.4321  5 

682.0858  7 

683.5654  7 

743.3526  6 

743.5596  6 

911.8574  8 

913.8633  8 

967.0898  7 

967.6143  7 

1172.680  9 

1175.290  9 

1223.021  8 

1223.940  8 

1464.327  10 

1467.614  10 

1510.683  9 

1512.071  5 

1786.498,  .  11 

1790.532  11 

1829.568  10 

1831.498  10 

2179.160  11 

2181.700  11 
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TABLE  XI 

Estimate  for  Fig.  4-5 
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/I  dz 

A(Jl  dz) 

^rot 

ATroc  ] 

CASE  1: 

97.5 

0.302 

155.3 

0.191 

CASE  2: 

115.1 

0.429 

156.1 

i 

0.232  j 

i 

CASE  3: 

115.1 

0.408 

158.1 

; 

0.225 

zt  -  88km  for  all  cases: 


TO 

atl 

AT0 

CASE  1: 

5.5 

5.5 

148.0 

0.3 

0.3 

CASE  2: 

6.0 

7.0 

148.0 

0.3 

0.3 

CASE  3: 

6.0 

7.0 

150.0 

0.3 

0.3 
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is  natural  because  the  least-square  fit  assumes  the  distribution 
similar  to  Gaussian. 

Geometric  Effect  for  Vertical  Observations 

The  discussion  so  far  devoted  to  the  "stationary  geometric 
effect,"  that  is,  the  geometric  effect  with  the  finite  thickness  layer 
at  rest.  Another  aspect  of  the  geometric  effect  is  the  "dynamic 
geometric  effect,"  becoming  the  basis  for  the  off-zenith  geometric 
effect  as  well  as  the  vertical  observation  with  the  emitting  layer 
moving  according  to  the  incoming  internal  gravity  wave.  This  is 
indeed  a  restatement  of  the  discussion  carried  out  by  Hines  and  other 
people  (7),  [11],  but  should  clarify  the  concept  of  the  dynamic 

geometric  effect.  Of  particular  interest  to  the  geometric  effect  are 
called  "polarization  relations"  by  Hines  [13],  because  they  define 
phase  as  well  <3  amplitude  relations  between  changing  physical 
quantities.  According  to  Hines,  the  air  parcel  travels  at  the  wind 
velocity  of  [13]: 


U  -  uKxC2  (kz — i  )  (4.21) 

C2 

W  -  u(w2-k2C2)  (4.22) 

where  U  and  W  are  phasers  for  horizontal  and  vertical  speed  respec¬ 
tively.  From  these  equations  it  is  obvious  that  u,  the  horizontal 
wind  speed,  and  w,  the  vertical  wind  speed,  have  a  phase  shift.  This 
makes  the  trajectory  of  the  air  parcel  an  ellipse  rather  than  a  line. 
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Since  there  is  also  a  common  factor,  exp(z/2H),  not  explicitly  repre¬ 
sented  in  the  polarization  relations,  the  true  trajectory  is  an  oval, 
asymmetric  ellipse  with  a  larger  radius  at  higher  altitude.  It  is 
this  unusual  air  parcel  trajectory  that  causes  the  dynamic  geometric 
effect,  both  in  vertical  and  off-zenith  observations.  The  geometric 
effect  for  8o-0  is  already  given  in  Chapter  II  and  III. 

Geometric  Effect  For  Off-Zenith  Observations 

The  extension  of  the  dynamic  models  introduced  so  far  allows  the 
simulation  of  the  observable  of  interest  when  the  line  of  sight  is  off 
to  the  zenith.  It  is  shown  that  a  small  modification  of  those  model 
is  enough  for  the  extension,  accounting  for  the  geometric  effect. 

There  is  a  simple  relationship  between  x  and  z  for  the  line  of 
sight.  When  the  line  of  sight  has  a  zenith  angle  80 ,  where  the  small 
o  indicates  'observation',  the  relation  between  x  and  z  along  the  line 
is : 


x  -  z  tan  80  (4.23) 


By  substituting  (4.23)  to  (4.4),  the  phase  factor  becomes: 


exp  i(wt-kxx-kzz)  -  exp  i(ut-K0z)  (4.24) 


where 


K0  -  -kxsin  So  +  kzcos  SQ  (4.25) 
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This  indicates  that  the  effect  of  the  off-zenith  observation  is 
the  apparent  wavelength  change  described  by  (A. 25).  It  should  be 
noted,  however,  that  the  derivation  of  (A. 25)  assumes  that  the  line  of 
sight  is  still  in  the  x-z  plane;  in  other  words,  parallel  to  the  wave 
vector  in  x-y  plane.  Further  geometric  considerations  allows  cases  of 
the  line  of  sight  not  to  be  parallel  to  the  wave  vector.  If  the  line 
of  sight  has  an  azimuth  angle  of  <p0  with  respect  to  the  wave  vector, 
the  equation  (A. 25)  is  modified  to  (refer  to  Fig.  A-7): 

K0  -  -kxsin  80  cos  <p0  +  k2cos  $0  (A. 26) 

Modifying  the  apparent  wavenumber,  or  the  apparent  wavelength, 
allows  the  previous  dynamic  models  for  the  vertical  observations  to  be 
used  in  the  same  manner  in  off-zenith  observation  simulation.  The 
only  change  occurs  when  multiplying  the  time  varying  term  (A.2A),  Ko 
should  be  used  instead  of  kx  and  kz.  Note  that  the  kx  and  kz  must 
still  be  used  to  calculate  other  parts  of  the  integrand,  because  those 
parts  are  determined  by  the  physically  meaningful  quantities  kx  and 
kz,  instead  of  the  which  is  only  a  result  of  the  geometric  effect. 
Therefore  K0  should  be  included  only  in  the  oscillation  term,  (3. A)  or 
(A . 2A) . 

Using  the  relation: 

tan  ta  -  -jr—  (A. 27) 
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where  $9  is  the  angle  of  the  wave  propagation  direction  to  the  zenith, 
(4.25)  and  (4.26)  becomes,  respectively: 

K„  -  k.cos  (1-  )  <‘.28) 

K0  -  k, cos  *.  (1-  Mn-£k  CM  ^  )  (4.29) 

A  simple  but  interesting  property  of  this  modifications  was 
briefly  introduced  in  Chapter  II.  When  8o-6a ,  equation  (4.29)  leads 
to  an  interesting  conclusion  when  v0“  0,  which  is  (4.28).  As  is  seen 
from  (4.25)  and  (4.26),  Ko-0,  the  apparent  wavelength  along  the  line 
of  sight  becomes  infinity  (Fig.  4-8).  Since  the  line  of  sight  is 
parallel  to  the  wave  plane  on  which  all  parts  have  the  same  phase,  the 
oscillation  is  observed  in  its  full  magnitude.  This  means  that  the 
geometric  effect,  the  cancellation  of  the  different  phase,  does  not 
occur,  since  there  is  no  phase  change  along  the  line  of  sight.  On  the 
other  hand,  any  geometrical  configuration  in  which  the  wave  In  the 
layer  along  the  line  of  sight  exercises  an  integer  multiple  of  360 
degree  phase  change,  as  in  the  simple  slab-layer  model,  no  osc'Mation 
will  be  observed  when  there  is  actually  an  oscillation.  In  ca=e  of 
more  realistic  layer  geometry  where  the  intensity  within  the  layer  is 
not  uniform,  such  cancellation  will  not  be  perfect  and  some  oscilla¬ 
tion  will  be  observed.  However,  serious  attenuation  in  apparent 
oscillation  amplitude  is  expected.  Note  that  the  oscillation  energy 
corresponding  to  those  attenuated  amplitude  does  not  disappear,  but  is 
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"transferred"  into  the  apparent  "DC"  component  of  the  emission,  or  a 
constant  brightness  representing  the  average  brightness.  This  effect 
is  especially  notable  with  an  infrared  imager  with  a  relatively  wide 
field  of  view.  The  same  wave  packet  with  the  same  physical  character¬ 
istics,  such  as  wavelength,  oscillation  frequency  and  amplitude,  will 
produce  entirely  different  signals  depending  upon  the  geometrical 
relation  between  the  line  of  sight  and  the  wave  propagation  direction. 
This  modification  of  the  true  wave  into  an  apparent  wave  can  be 
described  according  to  the  following  equation: 


Iobs^)  "  /  [0Ho  (z)  ]  exp(i(wt-K0z))  dz 
0 


(4.30) 


where  'Re'  denotes  the  real  part,  [OH0(z)]  is  the  stationary  con¬ 
centration  of  hydroxyl  at  altitude  z,  and  K0  is  given  by  (4.26)  or 
(4.29). 

When  [0Ho(z)]  is  approximated  with  a  constant  C  between  z-zL  and 
z-zu ,  zL<zu ,  and  0  elsewhere,  equation  (4.30)  becomes: 

C 

Iobs^)  “  ~  t  sin(ut-K0zL)  -  sin(wt-K0zu)  ]  (4.31) 

Ko 

2C  Mzu-Zl)  ko 

-  —  cos(  -  )  sin(wt - (z0+zL))  (4.32) 

Ko  2  2 

This  result  is  similar  to  the  oscillation  part  of  equation  (2.36). 
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The  replacement  of  the  true  wavelength  by  the  apparent  wavelength 
along  the  line  of  sight  is  necessary  and  sufficient  for  the  des- 
cription  of  the  off -zenith  geometric  effect. 
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CHAPTER  V 

DATA  PROCESSING  TECHNIQUE 

Dynamic  airglow  signal  processing  systems  intended  Cor  a  detection 
of  internal  gravity  waves  are  presented.  Selection  of  the  analysis 
method  and  choice  of  parameters  for  the  selected  method  depends  upon 
the  a  priori  knowledge  of  the  input  signal,  rather  than  the  signal¬ 
processing  theory  alone. 

Data  Reduction  System 

This  section  describes  data  of  the  mesospheric  infrared  radiation 
measured  by  the  interferometer  spectrometer.  A  description  and  the 
theory  of  operation  of  the  instrument  can  be  found  in  a  literature 
[33],  The  method  used  to  calculate  the  relative  emission  intensity 
and  the  rotational  temperature  using  the  observational  data  is 
described  in  Chapter  IV  and  Appendix  B.  Most  attention  is  focused  in 
an  introduction  of  a  method  to  calculate  IGW  parameters  such  as  time- 
resolved  power  spectral  density  out  of  the  observational  data.  These 
parameters  are  used  in  the  discussion  of  the  IGW  dynamics  and  the 
chemical  and  thermodynamic  modeling  of  the  infrared  emission  source, 
the  mesospheric  hydroxyl  layer. 

Michelson  Interferometer  Spectrometer 

The  interferometer  used  in  the  so-called  Fourier  transform 
spectroscopy  is  an  instrument  based  upon  the  same  interference 
principle  used  in  the  historical  Michelson-Morley  experiment  in  the 
1880's.  The  theory  of  operation  is  ^ound  in  literature  (33),  [34]. 
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Data  Reduction 


After  the  spectrum  of  the  incoming  infrared  light  is  calculated, 
the  emission  intensity  and  the  rotational  temperature  for  a  specified 
emission  band  are  calculated  as  functions  of  time.  In  the  present 
case,  developed  by  Espy  using  a  least-square  fitting  of  synthesized 
spectral  lines  of  one  vibrational  band  to  the  measured  spectral  lines 
of  the  same  band  was  used  [35].  This  method  also  provides  the  estima¬ 
tion  of  emission  intensity  and  error  estimates  of  these  values.  The 
least-square  method  is  described  in  Chapter  IV  and  Appendix  B  [31], 
[35].  These  values  are  analyzed  as  time  series  with  the  spectral- 
analysis  techniques.  Here,  the  word  "spectrum"  refers  to  the  spectrum 
in  the  temporal  oscillation  of  the  emission  intensity  and  the 
rotational  temperature,  rather  than  to  the  optical  emission  spectrum 
as  referred  to  previously. 

Spectral  Analysis  Techniques 

Introduction 

Once  the  temperature  and  intensity  data  are  prepared,  the  power 
spectral  densities  of  both  rotational  temperature  and  emission 
intensity  fluctuations  can  be  calculated.  These  frequency  domain 
spectra  are  useful  in  studying  theoretical  characteristics  of  IGW's  in 
their  frequency  domain  representations. 

One  frequency  domain  analysis  techniques  is  the  adoption  of  a 
Fourier  transformation  and  its  application,  a  short-time  Fourier 
analysis  [36].  Another  is  the  maximum  entropy  method  (MEM)  [37]. 
This  section  devotes  to  the  two  techniques. 
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Fourier  Analysis 

There  are  several  methods  of  estimating  the  power  spectral 
density.  One  commonly  used  method  is  calculating  the  Fourier 
transform  of  the  original  signal.  More  precisely,  the  power  spectral 
density  of  the  signal  is  obtained  by  the  Fourier  transform  of  the 
signal's  auto  correlation  [38].  This  method  does  not  provide  a  well 
defined  spectrum  when  the  original  signal  contains  only  a  few  cycles 
of  the  oscillation.  Most  experimental  data  contain  no  more  than  3  to 
7  cycles:  most  observations  last  45  minutes  to  a  few  hours  at  a  time, 
and  the  IGU's  to  be  studied  have  periods  of  tens  of  minutes  to  several 
hours.  Under  such  conditions,  the  best  estimate  of  the  oscillation 
frequency  has  the  uncertainty  Af  of: 

Af  -  ^2  (radian/second)  (3.1) 
Tl 

where  TL  is  the  total  data  length.  The  given  condition  makes  the 
value  of  the  peak  frequency  comparable  to  the  half  width  of  the  main 
lobe.  When  data  have  a  DC  component,  or  a  constant,  this  causes  an 
interference  of  DC  and  the  signal  of  interest.  Such  problems  may  be 
avoided  by  carefully  removing  DC,  as  are  discussed  in  the  following 
sections.  Also,  the  Fourier  analysis  fails  to  provide  a  pedk  when 
data  contain  less  than  one  cycle  of  oscillation.  In  such  case,  even 
DC  removal  does  not  work  properly. 

Also,  the  Fourier  transform  technique  produces  the  spectral 
function  with  significant  sidelobes  [38].  The  sidelobes  of  nearby 
frequencies  interfere  with  each  other,  making  the  Fourier  signal 
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analysis  inaccurate.  As  equation  (5.1)  indicates,  this  is  a  serious 
problem  when  the  data  contain  only  one  or  two  cycles .  In  such  case , 
the  sidelobe  of  the  data  window,  the  DC  component,  can  interfere  with 
the  main  lobe  of  the  signal. 

Short-Time  Fourier  Analysis 

It  is  necessary  to  use  a  signal  processing  technique  that  can  cope 
with  limitations  such  as  short  data  duration.  Moreover,  the  expected 
signal  is  not  monochromatic  sinusoidal.  Some  observational  data 
actually  show  such  chirp  oscillations,  as  shown  in  Fig.  5-1.  This  is 
real  IRFWI  data  taken  in  the  summer  of  1986  in  Fairbanks,  Alaska. 
The  expected  signal  does  not  contain  an  oscillation  of  a  single 
frequency  (monochromatic  wave),  so  the  limited  signal  duration 
problem  cannot  be  solved  even  with  any  instrumental  or  experiment 
operation  improvement,  such  as  longer  observation  time.  This 
invalidates  a  use  of  the  simple  Fourier  analysis. 

With  such  a  limitation,  one  technique,  a  short-time  Fourier  analy¬ 
sis,  can  be  an  appropriate  method  to  estimate  the  IGW  signal  chara¬ 
cteristics  [36],  This  technique  provides  a  means  to  estimate  instan¬ 
taneous  frequencies  constantly  changing.  This  method  is  described  in 
applicable  literature  [36].  With  this  technique,  a  spectrograph  of 
the  signal  can  be  derived.  "Spectrograph"  was  originally  used  in 
speech  analysis  as  a  graphic  representation  of  the  short-time  Fourier 
spectrum  [36].  When  possibility  exists  for  a  frequency  variation  in 
the  signal,  the  use  of  the  short-time  Fourier  analysis  may  be 
considered  with  careful  choice  of  the  window  length  which  should  be 
longer  than  the  signal  frequency  (period)  of  interest. 
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This  technique  is  stable  compared  to  the  maximum  entropy  method, 
especially  in  presence  of  noise.  It  will  be  shown  that,  unlike  the 
short-time  Fourier  analysis  technique,  the  maximum  entropy  method  is 
not  stable  in  presence  of  a  noise. 

Maximum  Entropy  Method 

A  technique  called  maximum  entropy  method  (MEM)  was  developed 
within  last  23  years  and  described  often  in  related  literature  [37J. 
This  method  has  the  following  advantage: 

1.  There  i:;  no  need  for  a  window  function.  MEM  provides  very  sharp 
peak  with  no  sidelobes.  This  is  due  to  the  fact  that  in  MEM,  the 
spectral  lines  are  approximated  by  a  function  of  the  form 
l/(e2+(w--w0)2)  where  e  is  a  parameter  [37]. 

2.  Very  high  frequency  resolution  can  be  achieved  using  very  small 
amounts  of  data,  often  less  than  a  complete  cycle. 

However,  MEM  has  the  following  problems  [37]: 

1.  An  interference  can  occur  between  spectral  lines.  Unlike  the 
similar  interference  due  to  sidelobes  in  the  Fourier  transform 
technique,  such  interference  produces  a  peak  at  a  wrong  frequency 
(frequency  shift)  in  MEM.  As  described  later  in  this  chapter, 
this  problem  may  be  avoided  by  use  of  an  appropriate  filter. 

2.  MEM  assumes  sharp,  or  narrow,  spectral  lines.  In  case  of  a  broad 
band  signal,  such  as  a  chirp  signal  which  does  not  have  a  sharp 
peak,  MEM  can  produce  faulty  spectrum  with  several  sharp  peaks. 
This  problem  can  sometimes  be  avoided  if  enough  high  order  is 
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used  in  the  calculation  of  MEM  spectrum.  When  the  amount  of  data 
is  too  small,  this  remedy  can  not  be  applied. 

3.  The  magnitude  of  the  resultant  spectrum  does  not  necessarily 

represent  any  physical  quantity  of  the  signal.  Even  their 

relative  magnitude  in  the  resultant  spectrum  may  not  represent 
their  relative  physical  magnitude.  A  better  mathematical 
understanding  of  MEM  is  needed. 

4.  MEM  is  susceptible  to  noise.  A  weak  noise  can  cause  splitting  of 

the  spectral  line  [37].  A  strong  noise  can  cause  a  shift  in  the 
spectral  line.  The  severity  of  such  effects  depends  the 

phase  and  duration  of  the  signal  in  data,  which  further  complica¬ 
tes  the  problem.  It  is  known  that  this  problem  can  not  be  cured 
by  increasing  the  amount  of  data. 

Several  results  of  real  data  processing  are  presented  in  the  last 
section  of  this  chapter.  Results  of  the  two  methods,  sho^t-cime 
Fourier  analysis  and  MEM,  are  compared. 

Digital  Filter 

Introduction 

The  interferometer  produces  one  interferogram  per  single  scanning 
of  the  moving  mirror,  or  the  prism  in  case  of  IRFWI.  The  airglow 
emission  intensity  and  the  rotational  temperature  are  calculated  for 
each  scanning  frame.  Because  of  the  way  data  are  reduced,  the  input 
signal  to  the  analysis  program  is  discrete  by  nature.  Therefore,  the 
application  of  the  digital  filter  technique  is  a  natural  choice  to 
design  filters  to  remove  the  effect  of  a  noise  and  other  interfere- 
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nces.  A  filter  is  needed  to  limit  the  signal  frequency  higher  than 
the  corresponding  frequency  of  the  window  length.  This  is  needed  for 
(short-time)  Fourier  analysis  to  remove  the  interference  of  the 
window.  As  will  be  shown,  MEM  can  also  produce  better  results  if 
enough  low  frequency  components  are  removed.  The  high  frequency 
removal  has  some  effect  on  MEM,  but  usually  not  in  Fourier  analysis. 
The  highest  IGW  frequency  component  of  interest  is  limited  by  the 
Brunt-Vaisala  frequency,  so  to  assist  a  time  domain  analysis,  a  use  of 
a  low  pass  filter  may  be  effective.  The  FIR  filter  technique  is  most 
suitable  in  processing  this  type  of  data. 

HR  and  FIR  Filter 

The  conventional  analog  filters  transformed  into  discrete  signal 
forms  are  called  infinite  impulse  response  filters  or  IIR  filters 
[39],  because  the  IIR  filter  uses  its  output  signal  as  the  feedback 
signal  to  the  input,  resulting  in  an  infinite  length  response  to  an 
impulse  input.  Because  of  the  feedback,  IIR  filters  are  sometimes 
called  recursive  filters  (39).  This  type  of  filter  has  the  advantage 
of  simple  design:  usually  one  can  design  this  type  of  filter  with  a 
calculator.  It  has  several  disadvantages  which  cause  serious  problems 
to  the  IGW  signal  analysis.  One  problem  is  the  nonlinear  phase 
characteristics  of  the  filter,  especially  near  the  cut-off  frequency. 
This  is  a  serious  disadvantage,  since  the  phase  difference  between  the 
intensity  and  the  temperature  oscillations  is  one  of  the  quantities  of 
interest.  Another  problem  is  that  an  IIR  filter  can  be  unstable. 
Because  the  IIR  has  a  feed  back  loop,  a  significant  impulse  response 
can  last  for  a  long  period  of  time.  This  can  introduce  a  faulty 
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oscillation  to  the  signal. 

Another  type  of  filter,  called  finite  impulse  response  filters  or 
FIR  filter,  can  be  designed  with  linear  phase  shift  characteristics 
while  maintaining  the  desired  amplitude  property  [39] ,  although  it 
usually  requires  a  large  number  of  data  points  to  obtain  the  superior 
phase  characteristics.  This  causes  truncations  of  data  on  both  ends. 
Since  the  number  of  data  points  from  the  instrument  are  usually  very 
limited,  the  FIR  filters  are  useful  only  when  enough  amount  of 
continuous  oscillations  are  recorded.  If  such  problems  are  overcome, 
FIR  filters  are  superior  to  HR  filters,  effective  in  studying  IGW 
oscillations  with  periods  of  30  minutes  or  less.  For  signals  with 
longer  periods,  a  simple  DC  removal  sometimes  works  better. 

The  effect  of  a  high-pass  filter  and  a  band-pass  filter  on  the 
maximum  entropy  method  is  demonstrated  in  Fig.  5-2.  Fig.  5-2(a)  is 
one  of  the  IRFWI  data  taken  in  Alaska.  Fig.  5-2(b)  is  several  spectra 
obtained  by  the  maximum  entropy  method.  The  high-pass  filter  has  a 
cut-off  frequency  represented  in  period  of  30  minutes.  The  band-pass 
filter  has  the  same  low  cut-off  frequency  and  a  high  cut-off  frequency 
of  3.3  minutes.  The  resultant  spectra  are  quite  different,  both  in 
amplitude  and  location  (frequency)  from  that  of  the  raw  data. 

An  effect  of  a  low  pass  filter  is  demonstrated  in  Fig.  5-3.  Fig. 
5-3(a)  i hows  a  sample  data  taken  by  IRFWI.  Fig.  5-3(b)  shows  Its 
spectrum  derived  by  MEM.  The  low-pass  filter  has  a  cut-off  frequency, 
represented  in  the  period,  of  3.3  minutes.  As  shown,  the  resultant 
spectra  are  quite  different  with  and  without  the  filter,  as  well  as  DC 
removal,  indicating  the  instability  of  the  MEM. 
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Another  use  of  the  low  pass  filter  is  demonstrated  in  Fig.  5-4. 
This  effect  is  a  suppression  of  the  apparent  oscillation  at  the 
sampling  interval  period  due  to  the  sampling  of  data  with  error  and 
noise.  Since  IGW  has  longer  period  than  Brunt-Vaisala  period 
theoretically,  such  spurious  spectrum  may  be  removed.  If  the 
differential  of  the  signal  is  of  interest,  for  instance,  numerical 
differentiation  can  give  good  results  only  after  such  noise  is 
removed.  Fig.  5-4(a)  is  another  data  from  Alaska.  Fig.  5-4(b)  is  the 
output  of  a  low  pass  filter.  Fig.  5-4(c)  shows  the  results  of 
numerical  differentiation  taken  from  the  original  and  filtered  signal. 

Time  Domain  Analysis 

Several  limitations  make  the  frequency  domain  analysis  of  the  IG'.C 
data  with  properties  such  as  short  duration  and  frequency  chirp 
inaccurate  and  difficult  to  implement.  These  limitations  are: 

1.  Missing  data  due  to  instrument  calibration  or  some  ocher  reasons. 
Data  are  truncated  during  these  periods.  Frequency  domain 
analysis  results  in  broadened  main  lobe  in  Fourier  spectrum  due  to 
averaging.  The  short-time  Fourier  technique  gives  lower  than  true 
frequency  in  some  cases. 

2.  Usually  only  few  cycles  of  oscillations  are  included  in  one  set  of 
data  even  without  any  data  truncation.  The  oscillation  period  of 
interest  is  usually  15  to  60  minutes  while  the  continuous 
measurement  duration  at  a  time  is  few  hours.  This  causes  a  large 
uncertainty  in  the  frequency  estimate. 

3.  Frequency  domain  analysis  does  not  always  properly  process  signals 
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86223:  Differentiation  without  filter 
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86223:  Differentiation  with  LPF  -3.3min 
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Fig.  5-4  Time  domain  analysis:  (c)  numerical  differentiation 
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observation  time,  do  not  occur,  because  the  original  wave  lasts 
only  a  few  periods  and  varies  its  ^frequency  rather  than  monochro¬ 
matic.  Use  of  the  short-time  Fourier  analysis  technique  offers  a 
limited  solution  to  these  problems.  MEM  may  provide  a  better 
solution  to  overcome  these  limitations,  but  its  own  problems  need 
to  be  resolved  before  MEM  can  be  applied  with  confidence. 

When  the  data  are  not  long  enough  and/or  signal  to  noise  ratio  is 
not  good  enough,  the  frequency  cannot  be  determined  from  the  data, 
even  if  the  signal  is  monochromatic.  Time  domain  analysis  would  work 
with  such  quality  of  data  with  fewer  problems  encountered  in  the 
frequency  domain  analysis.  Since  the  short-period  oscillation  can  be 
processed  with  the  short-time  Fourier  analysis  technique,  the  time 
domain  analysis  may  be  considered  to  process  oscillations  with  long 
period,  such  as  a  few  hours. 

In  the  time  domain  analysis  involving  long  period  oscillations, 
use  of  a  proper  low  pass  filter  is  essential.  This  is  especially 
true  with  the  rotational  temperature  data  with  smaller  signal  to 
noise  ratio  and  most  noise  has  high  frequency  components.  The  linear 
phase  shift  property  of  the  FIR  filter  is  very  significant  for  the 
time  domain  analysis  because  the  phase  relation  between  the  emission 
intensity  and  the  rotational  temperature  is  the  quantity  of  interest. 

The  output  of  a  low-pass  filter,  or  a  band-pass  filter  if  DC 
removal  is  needed,  can  be  examined  directly  in  time  domain.  Often  the 
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data  length  is  not  sufficient  to  deduce  the  frequency,  but  the 
existence  of  wave  activity  may  be  identified.  This  technique  may  be 
useful  when  there  is  a  large  and  slow  change  in  the  average  emission 
intensity  and/or  the  rotational  temperature  of  the  airglow. 

In  Fig.  5-4,  the  output  of  a  low-pass  filter  can  be  used  to 
determine  the  average  quantity.  Average  quantities  may  be  used  to 
evaluate  Krasovsky's  r)  or  the  background  stationary  quantity  in  a 
study  of  oscillations  with  much  shorter  periods. 

Example  Data  Analysis 

In  this  section,  the  short-time  Fourier  analysis  technique,  the 
MEM  and  the  time  domain  analysis  technique  are  applied  to  experimental 
data  taken  by  IRFWI  in  Pokerflat,  Alaska  during  the  summer  of  1986. 
The  dates  are  presented  in  their  universal  time  (UT)  format:  the  last 
two  digits  of  the  year  followed  by  three  digits  representing  the  day 
in  the  year  starting  January  1st  being  001.  Time  and  date  represents 
values  at  the  longitude  0.  The  emission  intensity  is  the  band 
intensity  for  (3,1)  band,  and  value  is  the  direct  output  from  the 
IRFWI  data  processing  system  developed  by  Espy  [35].  The  rotational 
temperature  is  reduced  by  the  same  system. 

Fig.  5-5 (a)  presents  the  raw  data  and  the  same  data  with  the 
linear  interpolations  and  the  simple  DC  removals.  Fig.  5-5(b) 
presents  their  Fourier  spectrum.  Fig.  5-5(c)  is  the  results  of  MEM  on 
these  data.  The  effect  of  a  simple  DC  removal  is  demonstrated  in 
these  examples.  The  DC  is  removed  by  simply  subtracting  the  average 
value  of  all  data  after  they  are  interpolated  with  a  straight  line 
connecting  the  edge  data  of  the  truncated  intervals. 
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86226  l:BPF:Hamming  85min 


86226  I:BPF:  MEM 


Fig.  5-5.  Signal  processing:  (f)  spectrum  of  output  of  band-pass 
filter  (BPF) 


146 


147 


The  spectrograph  can  be  plotted  from  any  direction.  The  viewing 
angle  of  the  plots  can  be  chosen  to  optimize  the  sight  for  spectral 
peaks  which  have  interesting  features.  Fig.  5-5(b)  contains  cne  such 
example,  indicated  by  an  arrow. 

Fig.  5-5 (d)  demonstrates  the  effect  of  a  high  pass  filter.  The 
comparison  reveals  a  major  difference  from  the  simple  DC  removal, 
especially  in  MEM.  This  figure  includes  the  combination  of  the  use  of 
the  filter  and  the  DC  removal,  in  this  order.  The  effect  of  the  DC 
removal  is  again  significant  in  MEM. 

Fig.  5-5 (e)  shows  the  effect  of  a  low-pass  filter.  Since  the  raw 
signal  contains  many  high  frequency  components,  the  effect  is  very 
significant  in  MEM.  The  DC  removal  causes  great  differences  in  all 
cases,  especially  with  MEM. 

Fig.  5- 5(f)  is  the  result  of  the  use  of  a  band  pass  filter.  The 
results  are  combinations  of  the  high  and  low  pass  filters.  Again,  the 
MEM  demonstrates  the  most  significant  effect. 

From  these  examples,  the  following  conclusions  about  the  signal 
processing  techniques  are  derived: 

1.  DC  removal  is  effective  in  all  cases  to  distinguish  long  period 
signals  and  minimizes  the  interference  of  the  DC  sidelobe  to 
short  period  signals. 

2.  Use  of  MEM  is  not  recommended  because  it  is  too  susceptible  to 
data  modification  such  as  filtering  and  DC  removal. 

3.  Long  period  signals  become  hard  to  detect  once  the  window  becomes 
shorter  than  the  signal  period.  Therefore,  window  length  should 
be  as  long  as  the  amount  of  data  allows. 
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4.  Hamming  windows  perform  reasonably  in  most  cases.  Use  of  Blackman 
windows  may  be  considered  when  the  interference  becomes  as  large 
as  the  signal  of  interest. 
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CHAPTER  VI 

DISCUSSION 


In  this  chapter,  results  of  both  theoretical  models  and  experimen¬ 
tal  techniques  are  analyzed  and  interpreted  .  First,  the  data 
processing  system  developed  in  Chapter  V  is  evaluated,  comparing 
performances  of  different  data  processing  techniques.  Interpretation 
of  the  chirp  signal  detected  by  the  signal  processing  system  using  the 
dynamic  IGW  model  is  demonstrated,  providing  the  basis  for  the 
conclusions  presented  in  Chapter  VII,  Secondly,  observables  for 
ground-based  optical  measurements  are  further  analyzed.  Finally,  the 
validity  and  utility  of  Krassovsky's  eta  parameter  are  discussed 
further,  and  the  extension  and  limitations  of  eta  are  discussed. 

Data  Processing  and  Interpretation 

In  Chapter  V,  several  signal  processing  techniques  were  intro¬ 
duced  and  applied  to  different  frequency  ranges.  These  techniques 
were  categorized  as  1)  high  (30  minutes  or  shorter),  2)  low  (30 
minutes  to  the  inverse  of  the  data  duration)  and  3)  trends  longer  than 
the  observation  window.  This  incrementation  is  based  upon  the 
relationships  between  the  oscillation  period,  the  data  length,  and  the 
window  length  used  in  the  short-time  Fourier  analysis.  These 
relationships  affect  both  the.  choice  of  an  appropriate  filter  and  the 
performance  of  the  signal  detection  technique. 

Accordingly,  this  section  is  divided  into  three  parts.  The  first 
part  treats  high  frequency  signal  processing;  the  second  part  is  the 
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low  frequency  signal  case,  and  the  third  part  is  the  long  trend  case. 
Design  specifications  of  the  finite  impulse  response  (FIR)  filters 
used  and  performances  of  the  short-time  Fourier  analysis  technique  and 
the  maximum  entropy  method  (MEM)  are  included.  Note  that  the  value 
of  30  minutes  used  to  differentiate  the  high  and  low  frequencies  is 
arbitrary  and  may  be  adjusted  according  to  the  data  length.  Longer 
data  streams  make  longer  boundaries,  45  minutes  or  more,  possible. 
This  boundary  value  is  chosen  to  distinguish  the  frequency  range  which 
exhibits  different  characteristics  of  the  interference  with  the  data 
window. 

High  Frequency  Signal  Case 

Relatively  high  frequency  signals,  with  a  period  of  30  minutes  or 
shorter,  can  be  detected  relatively  easily  compared  to  lower  frequency 
signals  because  the  data  duration  is  usually  much  longer  than  the 
signal  oscillation  period  so  that  a  suitable  window  length  for  the 
short-time  Fourier  analysis  can  be  easily  chosen.  The  same  condition 
also  makes  possible  straightforward  design  of  appropriate  filters, 
although  it  has  been  found  that  no  filter  is  necessary  in  most  cases 
to  detect  high  frequency  signal. 

Filters  are  not  needed  as  long  as  the  magnitude  of  the  signal  of 
interest  is  much  larger,  such  as  a  factor  of  10,  than  the  sidelobes  of 
the  undesired  frequency  peak  in  the  low  frequency  range.  When  a  low 
frequency  or  DC  signal  is  so  strong  that  its  sidelobes  may  interfere 
with  the  high  frequency  signal,  a  high  pass  filter  must  be  used.  Any 
filter,  including  an  FIR  filter,  has  some  ripples  in  its  pass  band 
characteristics,  so  the  use  of  a  filter  can  attenuate  the  desired 
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signal  as  well  as  the  unwanted  low  frequency  signal.  Results  with 
and  without  a  filter  need  to  be  compared  when  the  frequency  of  the 
signal  of  interest  is  close  to  the  filter  cut-off  frequency.  Also,  by 
using  windows  with  different  lengths,  which  produces  sidelobes  at 
different  frequencies,  the  distinction  between  a  sidelobe  and  a  faint 
signal  may  be  made.  It  is  within  this  frequency  range  that  most  of 
the  chirp  is  detected  with  the  short-time  Fourier  analysis  technique 
of  Chapter  V. 

In  case  of  high-pass  filters,  an  additional  DC  blocker  may  intro¬ 
duce  a  spurious  low  frequency  component  as  well  as  removing  DC,  but  a 
simple  clipping  of  such  large  components  revealed  that  they  do  not 
spoil  the  high  frequency  signal,  as  long  as  their  sidelobes  are  at 
least  10  times  smaller  than  the  signal  of  interest.  Since  an  I  GW 
spectrum  is  limited  by  the  Brunt-Vaisala  frequency,  an  appropriate 
band-pass  filter  may  be  used  for  visual  clarification  in  time  domain 
analysis  without  masking  possible  signals  produced  by  the  IGW.  The 
Brunt-Vaisala  frequency  in  the  mesosphere  corresponds  to  a  4  to  8 
minutes  period  depending  upon  the  temperature  (Appendix  A) . 

Although  use  of  a  filter  does  not  affect  the  output  of  the  short- 
time  Fourier  analysis,  it  alters  the  magnitude  of  spectral  peaks  for 
all  frequencies  in  MEM  output.  Also  in  MEM,  locations  of  some  peaks 
shift  as  a  result  of  filtering.  Because  of  such  instabilities,  the 
MEM  should  be  used  only  as  an  auxiliary  method  to  other  technique  such 
as  the  short-time  Fourier  technique.  Similar  frequency  shifts  and 
magnitude  variations  are  reported  as  a  result  of  noise  [37].  In  some 
case,  a  single  peak  splits  into  two  in  the  presence  of  weak  noise  for 
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reasons  not  fully  understood  [37], 

It  may  be  necessary  to  remove  the  DC  component,  slow  trends  and 
low  frequency  signals  when  these  signals  are  large,  causing  sig¬ 
nificant  sidelobes  in  high  frequency  region  which  can  directly 
interfere  with  high  frequency  signals.  With  the  use  of  an  appropriate 
window  function  which  produces  small  sidelobes,  this  difficulty  may  be 
avoided  without  removing  a  low  frequency  signal.  Fig.  6-l(a)  shows 
several  window  sidelobes  in  a  three-dimensional  spectrogram  of  the 
generated  chirp  signal  used  in  Chapter  III  (Fig.  3-8).  The  peak  of 
the  signal  has  a  magnitude  of  4,  and  it  is  noted  that  the  figure  has  a 
scale  of  order  of  10* 3 .  Comparison  with  the  example  in  Fig.  5-5 
suggests  that  the  small  peaks  in  Fig.  5-5(b)  are  real  signals,  not 
sidelobes  of  lower  frequency  signals.  These  peaks  are  almost  100 
times  larger  than  the  expected  sidelobe  magnitudes  for  the  Hamming 
window.  Without  such  a  window  function,  the  sidelobes  would  be 
comparable,  or  even  larger,  than  the  small  high-frequency  signals. 

Another  way  to  distinguish  the  sidelobe  of  a  low-frequency  signal 
and  a  true  signal  in  high-fr .quency  domain  is  the  use  of  windows  with 
different  length.  Figure  6-l(b)  represents  the  same  spectrograph  in 
Figure  5-5(b)  with  a  different  scale  of  signal  amplitude.  They  are 
the  same  data  with  three  different  windows,  60,  75  and  85  minutes 
long.  All  of  the  resultant  spectrograph  shows  small  peaks  in  UT-10-11 
at  usO.Ol.  The  magnitude  and  height  of  the  peaks  change,  but  the 
frequency  remains  the  same.  This  indicates  that  these  peaks  represent 
true  signal,  because  the  sidelobes  of  different'  length  windows  have 
different  frequencies. 
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Because  it  is  possible  to  detect  high  frequency  signals  under 
presence  of  rather  strong  low-frequency  signals  using  the  short- time 
Fourier  analysis  technique,  use  of  a  filter  may  be  avoided  since  a 
filter  truncates  the  data.  The  high  pass  filter  used  in  Fig.  5-5 
requires  180  data  points  to  calculate  an  output  value  and  truncates 
half  of  the  required  points,  90  points,  or  30  minutes  worth  of  data  at 
the  beginning  and  end  (20  second  sampling  interval  multiplied  by  90 
equals  30  minutes).  In  the  intensity  data  of  UT  date  86226,  the  use 
of  the  filter  causes  the  loss  of  signal  after  12  UT  which  contain 
significant  peaks  found  in  the  non- filtered  data. 

In  summary,  the  short-time  Fourier  analysis  technique  works 
without  any  major  problems  in  the  study  of  high  frequency  signals. 
The  FIR  filter  may  be  useful  in  time  domain  studies  but  is  not  needed 
in  short-time  Fourier  analysis.  Because  the  use  cf  filter  truncates 
both  ends  of  time  domain  data,  the  use  of  filter  should  be  avoided  in 
the  short-time  Fourier  analysis.  Chirp  signals  are  expected  to  change 
their  amplitude  as  well  as  frequency  as  functions  of  time,  so  careful 
study  of  the  spectrograph  is  needed  to  detect  them. 

Low  Frequency  Signal  Case 

Oscillations  with  longer  than  30  minute  periods  are  subject  to  a 
serious  interference  with  the  data  window.  Thus,  the  removal  of  DC  is 
necessary  and  can  be  accomplished  by  subtracting  the  average  from 
each  data.  MEM  sometimes  detects  a  peak  at  very  low  frequency,  even 
longer  than  the  data  period.  It  is  known  that.  MEM  is  capable  of 
detecting  signals  with  periods  longer  than  the  duration  of  the 
observed  data.  However,  in  the  presence  of  noise,  operations  such  as 
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the  removal  of  DC  or  use  of  filters  of  all  kinds ,  may  cause  faulty 
peaks  at  all  frequencies  (37],  and  therefore  may  not  be  completely 
accurate.  The  use  of  the  time  domain  analysis  may  be  used  to  test  the 
result  of  MEM. 

Ir  more  than  one  signal  cycle  is  included  in  the  data,  with  DC 
removal,  the  short-time  Fourier  analysis  technique  can  be  used.  Since 
the  uncertainty  in  frequency  is  comparable  with  the  oscillation 
frequency,  the  detection  of  the  frequency  variation  as  a  result  of 
chirp  may  not  be  performed  well  in  this  frequency  range. 

As  the  example  in  Fig.  6-l(a)  indicates,  the  resultant  spectri.  of 
the  short-time  Fourier  analysis  varies  according  to  the  phase  of  the 
signal  within  the  window.  Similar  effects  of  the  signal  phase  is 
reported  for  MEM  [37].  The  time  domain  signal  analysis  is  needed  to 
assess  the  result  of  frequency  domain  analysis  by  observing  such  phase 
characteristics . 

Trends 

In  this  study,  the  word  "trend"  refers  to  a  slow  signal  variation 
caused  by  a  very  long  period  oscillation  or  a  slow  change  due  to  some 
other  factor.  The  Fourier  transform  technique  may  not  be  applied  to  a 
signal  when  data  contain  less  than  a  single  cycle  of  the  signal.  MEM 
is  capable  of  detecting  a  peak  in  such  cases,  but  the  resultant 
spectra  are  subject  to  noise  and  the  phase  of  signal  in  data  [37]. 
The  effect  of  filtering  may  be  interpreted  as  a  result  of  redistribu¬ 
tion  of  noise  by  such  filters.  Therefore,  the  resylts  of  MEM  analysis 
need  to  be  analyzed  with  an  aid  of  some  other  method.  One  such 
method  used  in  this  dissertation  is  the  representation  of  signal  in 
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the  time  domain.  The  low  pass  filter  may  be  used  to  clarify  the  low 
frequency  variation  (Fig.  5-4).  Fig.  6-2  represents  the  MEM  spectrum 
of  the  same  data,  detecting  a  peak  at  about  70  minutes.  The  visible 
trend  in  the  time  domain  representation  is  slightly  longer  than  an 
hour,  giving  consistent  results  (Fig.  5-4). 

It  is  noted  that  whether  a  signal  frequency  is  "low"  or  "high"  is 
determined  only  with  respect  to  the  window  length,  and  not  related  to 
any  physical  characteristics  of  IGWs.  Several  windows  with  various 
lengths  may  be  used  for  different  frequency  range  analysis.  The  three 
dimensional  spectrograph  will  help  such  overall  approach. 

Source  Distance  Estimation 

The  distance  between  the  source  of  an  observed  IGW  and  the 
observed  location  can  be  estimated  using  the  spectrograph  and  the 
dynamic  IGW  model  presented  in  Chapter  III. 

Fig.  6-3  shows  two  days  of  IRFWI  data  taken  in  Pokerflat,  Alaska 
(147.5  W  65.1  N) ,  during  the  summer  of  1986.  Fig.  6-3(a)  and  (b)  are 
spectra  of  the  intensity  and  the  temperature  of  the  UT  day  216,  and 
Fig.  S-3(c)  and  (d)  are  spectra  of  the  intensity  and  the  temperature 
of  the  UT  day  '.23  respectively. 

These  examples  contain  several  peaks  indicated  by  the  arrows  in 
which  frequencies  change  with  time.  It  is  expected  in  the  result 
obtained  in  Chapter  III  that  the  high  frequency  IGW  si^-ials  are  likely 
to  chirp  towards  higher  frequency.  The  series  of  peaks  in  the 
spectrograph  of  the  temperature  data  on  day  223  are  consistent  wicn 
the  model  prediction.  Signals  longer  than  one  hour  are  likely  to 
chirp  towards  a  lower  frequency  at  a  rate  indicated  in  Fig.  3-11. 
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Fig.  6-3.  Fourier  spectrum:  (a)  Intensity:  day  86216 
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Fig.  6-3.  Fourier  spectrum:  (b)  Temperature:  day  86216 
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Fig.  6>3.  Fourier  spectrum:  <c)  Intensity:  day  86223 
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Fig,  6-3.  Fourier  spectrum:  (d)  Temperature:  day  86223 
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The  interference  of  the  window  and  low  frequency  signal,  as  well  as 
such  frequency  change,  makes  the  detection  of  a  chirp  at  low  frequency 
difficult. 

The  temperature  contains  a  more  complex  spectral  structure  than 
the  intensity  due  to  the  geometric  effect  explained  in  Chapter  IV, 
where  it  was  demonstrated  that  the  observable  intensity  is  subject  to 
different  cancellation  mechanisms  from  the  observable  rotational 
temperature,  resulting  in  more  serious  cancellations.  In  the  day  223 
measurement,  the  large,  low  frequency  peak  in  the  intensity  spectrum 
at  9.6  UT  is  not  observed  in  the  temperature  spectrum.  This  can  be 
also  seen  in  the  time  domain  representation  in  Fig.  6-2(a).  The 
temperature  curve  is  almost  flat  after  rising  significantly.  This  may 
not  be  a  part  of  a  sinusoidal  wave. 

According  to  Fig.  3-11,  for  each  IGW  elevation  angle  a,  certain 
frequency  components  travel  at  the  fastest  speed,  and  the  chirp  starts 
towards  both  high  and  low  frequency  from  this  particular  frequency. 
The  low  frequency  component  has  an  asymptotic  value  given  by  equation 
(3.7)  and  Fig.  3-9  for  each  a.  The  high  frequency  limit  is  always 
the  Brunt-Vaisala  frequency.  Near  the  fastest  travelling  frequency, 
the  frequency  variation  rate  is  much  larger  than  near  the  both 
limiting  frequency  (Fig.  3-11). 

It  is  possible  to  estimate  a,  therefore  the  distance  of  the  IGW 
source,  using  the  frequency  variation  rate.  A  fast  frequency 
variation  indicates  that  the  observed  frequency  is  close  in  value  to 
the  fastest  travelling  frequency.  A  slow  frequency  change  rate 
observed  in  the  low  frequency  branch  indicates  that  the  observed 
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frequency  is  close  in  value  to  the  limiting  frequency.  By  estimating 
either  the  fastest  travelling  frequency  or  the  limiting  frequency,  the 
elevation  angle  a  may  be  estimated.  This  is  a  useful  technique  by 
which  the  distance  of  the  IGW  source  may  be  estimated.  This  model  may 
be  tested  by  measuring  both  the  fastest  and  limiting  frequencies  and 
testing  if  the  a  obtained  by  the  two  methods  are  consistent. 

There  is  another  way  to  estimate  a  if  a  chirp  signal  is  found  in 
the  frequency  range  between  both  limiting  frequencies  and  near  the 
fastest  travelling  frequency.  The  fastest  travelling  frequency  can  be 
determined  whether  it  is  lower  or  higher  than  the  observed  frequencies 
by  the  direction  of  the  frequency  chirp,  as  shown  in  Fig.  3-11.  This 
also  allows  for  an  estimate  of  the  fastest  travelling  frequency. 
Therefore,  the  distance  of  the  IGW  source  can  be  estimated. 

In  the  example  of  Fig.  6-3(a)  and  (b),  the  day  86216,  the  three 
peaks  at  UT-9  to  10  and  w=0.C26,  0.016  and  0.012  (rad/sec),  or  T=4, 
6.5  and  8.5  (minutes),  exhibit  chirp  towards  higher  frequency.  This 
example  also  has  a  peak  at  l'T-8.8  to  9.5  and  has  w* 0.006,  or  T=17 
(minutes),  chirping  towards  a  lower  frequency  which  is  visible  in  the 
temperature  spectrograph  but  net  in  the  intensity.  According  to  Table 
VII  and  Fig.  3-9,  this  indicates  that  a  is  less  than  20  degrees,  and 
the  distance  to  the  source  can  be  estimated  longer  than  240  kilo¬ 
meters  . 

In  the  example  of  Fig.  6-3(c)  and  (d) ,  on  the  day  86223,  a  few 
small  peaks  at  high  frequency  are  also  seen,  but  there  are  two  large 
peaks  of  particular  interest  (only  one  is  visible  in  the  intensity 
spectrograph).  They  appear  as  the  fastest  travelling  frequency 
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components  followed  by  chirps  towards  both  high  and  low  frequencies. 
Their  fastest  travelling  frequency  components  are  0.008  and  0.005 
(rad/sec),  or  T-13  and  21  minutes.  They  correspond  to,  according  to 
equation  (3.7)  and  Table  VII,  a-28  and  17  (deg)  and  1^-170  and  290 
(km) . 

Such  distance  estimation  will  help  to  identify  the  source, 
especially  if  combined  with  the  directional  information.  In  the  above 
examples,  the  direction  was  not  available. 

Observables 

The  geometric  effect  significantly  influences  observed  emission 
oscillations,  especially  in  the  high  frequency  components  which  is 
associated  with  short  vertical  wavelengths.  The  effect  on  the 
observed  rotational  temperature  can  be  significant,  but  the  magnitude 
of  the  effect  was  found  to  be  on  the  same  order  of  uncertainties 
caused  by  other  mechanisms  such  as  instrument  error,  noise  in 
instrument,  noise  in  signal  and  water  vapor  absorption.  The 

geometric  effect  upon  the  rotational  temperature  is  two-fold:  1) 

cancellation  in  oscillation  amplitude  and  2)  non-gaussian  distribu¬ 
tion  of  brightness -weighted  temperature.  Further  study  is  needed  to 
reduce,  or  if  possible  to  remove,  these  uncertainties. 

Using  the  off-zenith  geometric  effect  model,  it  is  possible  to 
estimate  IGW  parameters  with  a  measurement  covering  a  wide  viewing 
angle.  The  IGW  parameters  such  as  u,  kx  and  kz  can  be  reduced  using 
the  off-zenith  geometric  effect  model,  provided,  that  the  apparent 
wavelength  can  be  measured.  To  observe  the  wavelength,  it  is 
necessary  to  measure  a  wide  field  of  view  (fov)  that  extends  more 
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than  one  apparent  spatial  wavelength  at  sufficient  resolution.  The 
resolution  must  be  at  least  half  the  apparent  wavelength  to  satisfy 
the  sampling  theorem.  Because  the  direction  of  propagation  needs  to 
be  measured  at  the  same  time,  to  avoid  aliasing,  the  resolution  within 
the  fov  should  be  much  higher  than  the  Nyquist  frequency.  For  an  IGW 
with  the  period  of  15  minutes  to  several  hours,  the  theoretical 
apparent  wavelength  is  between  a  few  tens  to  a  few  hundred  kilometers. 
The  packet  speed  of  such  wave  is  150  to  250  m/s,  using  equations 
(A. 14)  and  (A. 15),  so  it  takes  about  one  hour  for  the  wave  to  travel 
on  the  order  of  500  kilometers.  It  is  necessary  to  measure  more  than 
several  periods  of  a  wave  without  major  interruptions.  One  measure¬ 
ment  instrument  which  can  satisfy  such  requirements  is  an  infrared 
imager  [5],  The  use  of  an  infrared  all-sky  camera  will  help  evaluate 
data  from  instruments  with  smaller  fov.  As  the  viewing  elevation 
angle  decreases,  closer  to  the  horizon,  the  optical  path  becomes  more 
involved  with  air  refraction,  scattering  and  water  vapor  absorption, 
as  well  as  more  serious  geometric  effect.  Therefore,  the  low-angle 
observation  should  be  limited  to  the  emission  intensity  measurement, 
not  the  rotational  temperature  which  is  more  sensitive  to  those 
interferences . 

Using  the  Hines'  IGW  model,  a  relation  between  w,  kx ,  kz  and  the 
distance  from  the  IGW  source  and  the  observation  location  can  be 
obtained.  Equation  (3.6)  gives  the  kx ,  or  the  horizontal  wave  length 
as  a  function  of  w  and  the  elevation  angle  a  of  the  wave  propagation. 
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U4  (u2  -  Xw|  ) 

k2  -  —  -  (3.6) 

C2  (u2-w2)(o2  +  (w2-<j2)tan2  a) 

Assuming  Che  heighc  of  che  airglow  layer,  Che  distance  of  the  wave 
source  can  be  calculated  using  the  elevation  angle.  Fig.  6-4  gives 
such  relations  for  several  oscillation  period.  Since  the  observed 
vertical  wavelength  is  usually  8  to  15  km  regardless  of  frequency,  it 
seems  that  relatively  small  portion  of  possible  waves  are  actually 
observed.  From  this  relation,  it  may  be  possible  to  estimate  the 
distance  of  the  source  using  measured  u  and  kx  or  kt .  Because  Xx  ,  the 
horizontal  wavelength,  is  a  steep  function  of  the  distance  and  u,  it 
is  practical  to  test  the  correspondence  between  w  and  Ax  first. 
Because  the  observed  vertical  wavelength  is  limited,  the  distance  may¬ 
be  obtained  from  either  w  or  Ax  alone  since  the  two  have  almost  one  to 
one  relation.  More  observational  data  including  wavelengths  need  to 
be  accumulated  to  test  this  method. 

Discussion  on  Eta 

In  this  section,  the  implementation  and  utility  of  Krassovsky's 
parameter  rj  is  further  discussed.  As  was  shown  in  previous  chapters, 
experimental  data  would  seldom,  if  at  all,  be  appropriate  to  calculate 
eta  as  originally  defined.  Instead,  signal  processing  techniques 
suitable  for  the  likely  data  characteristics  were  developed. 

The  original  t)  is  defined  under  the  assumption  that  the  airglow 
layer  is  much  thinner  than  the  vertical  wavelength  of  an  IGW. 
Therefore,  in  case  of  long  period  oscillations  with  long  wavelengths, 
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such  as  a  few  hours  with  longer  than  500km  horizontal  wavelength,  the 
rj  may  still  be  calculated  using  the  original  definition.  Such 
observational  data  conditions  require  several  hours  of  continuous 
measurement  with  the  use  of  a  low-pass  filter.  The  applicability  of 
the  time  domain  signal  processing  is  most  suitable  in  this  kind  of 
study,  since  it  is  a  similar  method  used  by  Krassovsky  in  his 
introduction  of  rj . 

The  concept  of  17  may  be  utilized  in  the  case  where  rj  is  not 
limited  to  a  quasi-sinusoidal  wave  with  constant  background  for  both 
the  emission  intensity  and  temperature  As  was  suggested  in  Chapter 
II,  I  and  T  can  be  instantaneous  values  and  dl  and  dT  can  be  the  time 
derivatives  of  the  emission  intensity  and  the  temperature  respective¬ 
ly.  The  resultant  rj  should  be  of  the  same  value  according  to  (2.53). 
The  derivation  by  Krassovsky  did  not  include  the  geometric  effect 
introduced  in  the  present  study,  and  as  a  consequence,  (2.53)  is 
limited  in  accuracy  of  values  of  tj  in  model  calculations. 

In  computing  both  intensity  and  temperature,  as  is  shown  in  Fig. 
6-5(a),  each  data  point  has  an  error,  or  uncertainty.  The  central 
values  exhibit  discontinuities,  resulting  in  the  large  signal 
component  at  the  frequency  corresponding  to  the  time  interval  between 
data  points.  Since  this  comp«  nt  caused  by  the  sampling  is  spurious, 
an  optimal  filtering  process  is  necessary.  Because  of  the  phase 
property,  an  FIR  digital  filt  :  is  desirable,  provided  that  sufficient 
data  p  >ints  are  available.  Without  filtering,  calculating  dT  as: 

dT  -  T ( t+A  t ) -T ( t )  (6.1) 
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will  result  in  a  large  value  for  dT  that  renders  such  calculations 
useless  for  evaluating  r) .  Some  of  the  outputs  of  low-pass  filters  and 
their  numerical  differentiations  are  given  in  Fig.  6-5(b). 

Another  consideration  is,  as  expected  from  the  theoretical 
calculation  and  some  observed  value  of  t) ,  that  the  fluctuations  in 
temperature,  dT,  are  usually  very  small.  The  observational  77  is 
reported  to  be  between  5  to  10  in  many  cases  [4],  [6].  This  makes 
dT/T  much  smaller  than  dl/I  so  that  expected  value  for  dT  is  often 
undetectably  small.  Quite  often,  possible  fluctuations  in  temperature 
are  smaller  than  the  error  estimate,  or  even  the  small  jump  in  each 
data  point  caused  by  discrete  measurements.  In  such  cases,  because 
the  signal  can  be  smaller  than  noise,  the  filter  eliminates  not  only 
noise  but  also  the  signal.  Because  the  frequency  components  of  the 
sampling  induced  noise  and  of  the  noise  associated  with  the  error 
estimate  are  on  the  same  order  as  the  sampling  frequency,  long  period 
fluctuations  could  be  detected  with  an  appropriate  low-pass  filter. 
In  case  of  a  sampling  period  of  about  20  seconds,  long-period 
fluctuations,  s"ch  as  hours,  could  be  detected,  provided  that  at  least 
one  entire  cycle  is  recorded. 

Most  observations  report  the  eta  to  be  about  5  to  10  [4],  [6]. 
Thus  the  intensity  fluctuations  are  much  larger  than  the  temperature 
fluctuations,  so  this  problem  is  at  least  5  to  10  times  less  serious 
in  the  intensity  data.  One  possible  method  to  overcome  this  problem 
in  temperature  data  is  using  the  filtered  intensity  signal  as  a 
reference  to  detect  temperature  fluctuations  using,  say,  the  coherent 
detection.  Since  a  theoretical  relation  between  observed  intensity 
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and  temperature  is  available,  comparison  between  expected  and  observed 
temperature  fluctuations  using  such  model  may  be  possible. 

As  the  short-time  Fourier  analysis  revealed,  the  temperature  data 
exhibits  more  oscillations  of  a  variety  of  frequencies  than  the 
intensity.  The  theoretical  model,  the  geometric  effect,  also  predicts 
that  the  temperature  and  the  intensity  are  subject  to  different 
cancellation  mechanisms.  It  was  shown  in  Chapter  II  that  when  there 
is  a  temperature  variation  within  the  emission  layer,  temperature  may 
exhibit  some  oscillation  while  the  intensity  oscillation  is  attenuated 
significantly.  Such  conditions  also  make  17  less  useful,  particularly 
for  IGW's  with  short  vertical  wavelengths,  namely,  comparable  with  the 
airglow  layer  thickness. 

Also,  not  all  fluctuations  are  necessarily  caused  by  IGWs.  The 
characteristics  of  the  airglow  fluctuations  caused  by  IGWs  can  be 
predicted  using  theoretical  models,  as  is  introduced  in  previous 
chapters.  Therefore,  fluctuations  in  observed  data  that  do  not 
satisfy  those  properties  could  be  considered  to  have  other  origin. 
The  long  period  trend  need  not  be  a  part  of  a  sinusoidal  oscillation, 
for  instance.  This  paper  does  not  deal  with  any  of  the  airglow 
fluctuations  caused  by  something  other  than  the  linear  IGW's. 
However,  the  result  of  this  study  can  be  used  to  reject  the  pos¬ 
sibility  of  IGW- induced  fluctuations  on  certain  data  which  do  not 
satisfy  the  predicted  characteristics. 
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CHAPTER  VII 


CONCLUSIONS 


Summary 

The  major  achievements  of  this  study  are  summarized  as  follows: 

1.  Developed  a  systematic  formulation  of  a  photochemical  model  of  the 
mesospheric  hydroxyl  layer.  The  model  can  be  used  to  calculate 
airglow  emission  line  intensities  with  IGW  parameter  inputs  for 
comparison  with  airglow  observations. 

2 .  Summarized  the  dynamics  of  the  IGW  in  a  form  appropriate  for  the 
estimation  of  ground-based  observations  of  mesospheric  airglow 
emissions.  It  was  found  with  a  realistic  model  atmosphere  that 
the  fluctuations  in  the  airglow  emissions  do  not  necessarily 
represent  the  amplitude  nor  the  phase  of  the  IGW  oscillation. 
Only  the  wave  frequency  is  represented  correctly. 

3.  Found  that  the  IGW  propagation  elevation  angle  is  limited  by  a 
function  of  the  IGW  frequency.  With  this  finding,  it  is  possible 
to  estimate  a  minimum  distance  of  the  source  of  an  observed  IGW 
using  its  frequency. 

4.  Found  that  IGW's  are  likely  to  be  quasi-sinusoidal ,  often  referred 
to  as  "chirp",  after  travelling  a  long  distance  into  the  meso¬ 
spheric  region.  As  a  result,  IGW- induced  fluctuations  in  the 
mesospheric  airglow  emissions  are  likely  to  be  variable  frequency 
sinusoidal.  This  is  similar  to  whistler  mode  electromagnetic  wave 
propagation  in  the  ionosphere.  However,  the  changing  rate  of  the 
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frequency  was  found  to  be  small,  so  observations  exhibiting  a 
cycle  or  two  oscillation  can  be  treated  as  monochromatically 
sinusoidal . 

5.  Summarized  the  geometric  effect  of  ground-based  optical  observa¬ 
tion..  This  effect  was  used  in  conjunction  with  the  dynamics  of 
the  IGW  to  evaluate  the  zenith-aspect  observational  data.  Tc  was 
found  that  the  brightness -weighted  kinetic  temperature,  integrated 
along  the  line  of  sight,  is  likely  to  be  depart  from  a  gaussian 
distribution.  Since  the  rotational  temperature  is  related  to  the 
kinetic  temperature,  this  departure  affects  the  apparent  rotation¬ 
al  temperature.  Use  of  the  least  square  fitting  technique  to 
estimate  the  rotational  temperature  is  therefore  in  question.  A 
systematic  error  in  the  rotational  temperature  estimation  is 
likely  to  occur;  however,  this  error  is  usually  less  than  1  K 
when  using  rotational  lines  with  the  rotational  level  j<6.  This 
would  be  significant  for  rotational  temperature  estimates  of  about 
1  K  or  less  uncertainty.  As  the  temperature  estimate  accuracy  and 
precision  are  improved,  the  geometric  effect  on  rotational  lines 
will  become  increasingly  significant. 

6.  Analyzed  off -zenith  observations  as  an  extension  of  the  zenith 
observation  geometric  effect.  The  result  can  be  used  to  optimize 
both  experimental  planning  and  data  analysis.  The  off-zenith 
observation  can  be  unambiguous  only  if  the  background  conditions, 
such  as  the  stationary  kinetic  temperature  and  the  concentrations 
of  related  chemical  species,  are  uniform  horizontally  along  the 
line  of  sight. 
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7.  Developed  signal  processing  techniques  optimized  for  the  detection 
of  IGW- induced  chirp  oscillation.  The  expected  chirp -like 
frequency  change  appears  to  occur  frequently,  and  the  short  time 
Fourier  transform  technique  works  without  problem  when  analyzing 
data  containing  more  than  two  complete  cycles  of  oscillation.  As 
the  amount  of  data  becomes  smaller,  approaching  only  a  single 
oscillation  period,  the  technique  gives  lower  than  actual 
frequency.  When  measurement  window  contains  less  than  one 
oscillation  period,  this  technique  does  not  detect  a  spectrum  peak 
for  that  signal. 

8.  Found  that  the  Krassovsky's  rj  is  useful  only  to  limited  cases. 
Such  casis  need  to  satisfy  the  following  conditions: 

a)  a  stationary  background  in  emission  intensity  and  temperature, 

b)  several  cycles  of  oscillations  are  recorded  in  both  emission 
intensity  and  rotational  temperature  data  and  c)  noise  and  error 
estimates  of  both  data  are  smaller  than  their  signal  amplitude. 

Conclusions 

The  following  conclusions  are  derived  from  this  study: 

1.  It  was  found  that  the  temperature  dependence  of  the  rates  of  the 
chemical  reactions  responsible  for  the  formation  of  the  hydroxyl 
layer  can  have  a  significant  effect  in  the  translation  of  the  IGW 
and  the  observable  in  the  passive  ground-based  optical  measure¬ 
ment.  The  significance  of  this  effect  varies,  and  can  only  be 
evaluated  with  knowledge  of  the  vertical  distribution  of  physical 
parameters,  such  as  concentrations  of  related  chemical  species  and 
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temperature.  As  a  result  of  this  effect,  oscillations  observed  in 
hydroxyl  airglow  do  not  accurately  represent  IGW  amplitude. 

2.  The  distance  between  the  source  and  the  observation  site  of  an 
IGW  may  be  deduced  from  the  observed  oscillation  frequency  and  its 
changing  rate.  The  estimation  becomes  more  accurate  if  the  low 
limiting  frequency  of  the  chirp  is  observed,  since  the  low 
limiting  frequency  is  a  unique  function  of  the  IGW  propagation 
elevation  angle. 

3.  The  relation  between  frequency  and  horizontal  wavelength  may  be 
used  to  estimate  the  distance  to  the  source. 

4.  If  the  physical  parameters  such  as  temperatures  along  the  IGW  path 
are  known,  the  spectrum  of  the  wave  source  may  be  estimated  using 
the  observed  oscillation  spectrum  and  wavelengths.  In  turn,  if 
the  spectrum  of  the  source  is  known,  such  physical  parameters 
along  the  wave  path  may  be  estimated  using  observational  data. 

5.  Off-zenith  observations  can  be  evaluated  using  a  geometric  effect 
model.  According  to  the  model,  aiming  at  a  particular  angle 
implicitly  sets  up  spatial  filtering  due  to  the  slanted  nature  of 
the  IGW  phase  plane.  Waves  are  attenuated  or  amplified  according 
to  their  frequency,  wavelength  and  propagation  direction  with 
respect  to  the  line  of  sight  for  each  observation  angle,  resulting 
in  favored  range  in  frequencies  and  wavelengths  for  that  direc¬ 
tion,  similar  to  a  band-pass  filter. 

6.  The  use  of  the  maximum  entropy  method  (MEM)  is  not  recommended  in 
a  study  of  the  mesospheric  airglow  oscillations. 

7.  The  application  of  Krassovsky's  ij  upon  experimental  data  has 
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several  limitations.  Such  limitations  are:  1)  the  -.umerical 

differentiation  of  data  is  heavily  affected  by  the  characteristics 
of  the  preconditioning  filters,  2)  the  temperature  fluctuation  is 
usually  very  small,  making  ij  very  sensitive  to  a  numerical  error, 
because  this  small  value  is  in  denominator. 

8.  The  application  of  Krassovsky’s  t]  may  be  appropriate  for  some 
cases.  Such  cases  include  long  period  IGWs,  more  than  several 
hours  to  a  few  days,  which  tend  to  have  vertical  wavelengths  much 
longer  than  the  hydroxyl  layer  thickness,  approaching  to  the 
assumption  in  the  definition  of  the  rj. 

9.  Off -zenith  angle  observation  may  help  in  some  cases  where  such 
geometrical  configuration  causes  amplifications  to  some  range  of 
signals.  This  gain  has  a  penalty  of  attenuations  to  signals 
outside  of  the  gain  range,  due  to  the  same  geometric  effect.  This 
range  is  a  function  of  the  viewing  elevation  angle . 

Recommendations 

Based  on  the  findings  of  this  study,  the  following  recommenda¬ 
tions  are  presented: 

1.  Because  the  volume  emission  rate  varies  over  the  hydroxyl  layer, 
future  studies  require  depth  data  to  be  obtained  in  order  to 
compare  the  theoretical  model  with  experimental  data.  Some 
methods  to  obtain  such  data  are:  a)  use  of  ranging  instruments, 
such  as  RADAR  (Radio  Detection  And  Ranging)  and  LIDAR  (Light 
Detection  And  Ranging),  b)  stereoscopic  observations  by  use  of  a 
single  instrument  which  scans  the  sky  or  multiple  count  instru- 
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merits  at  distant  locations.  The  recommendation  for  the  range  of 
each  instrument  is  presented  in  part  4.  The  distance  between 
different  instruments  for  multiple  site  observation  needs  to  be 
large  enough  for  the  geometric  effect  to  be  observable.  This  is 
related  to  the  range  of  each  instrument.  Separation  of  100  to 
1000  km  is  recommended  for  observation  of  oscillations  over  a  few 
hour  period.  These  waves  have  low  propagation  elevation  angle, 
10  degrees  or  fewer,  so  the  longer  distance  may  be  adopted. 
Shorter  period  oscillations  may  need  to  be  observed  with  shorter 
separation,  such  as  100  to  300  km,  because  such  wave  may  propagate 
at  larger  elevation  angle,  30  degrees  or  so.  With  the  acquisi¬ 
tion  of  stereoscopic  data,  the  wave  dynamic  model  can  be  tested 
more  accurately,  and  the  location  of  the  source  may  be  calculated. 

2.  Once  the  dynamic  model  is  verified,  some  depth  information  may  be 
deduced  using  a  simple  ground-based  passive  observation  data  and 
the  geometric  effect  model.  Therefore,  further  improvement  of  the 
I GW  dynamic  model,  especially  nonlinear  model,  is  recommended. 

3.  Some  of  the  techniques  for  modeling  the  geometric  effect  may  be 
applied  to  a  modeling  of  other  types  of  observations,  such  as 
satellite-based  passive  optical  observation.  Further  study  in 
this  area  is  recommended. 

4.  Because  IGWs  are  likely  to  be  in  the  form  of  spatially  limited 
chirp  wave  packets,  observation  of  the  same  wave  packet  as  it 
propagates  through  the  atmosphere  for  long  period  of  time  is 
preferred  over  the  observation  of  a  fixed  direction  which  provides 
data  on  a  small  fraction  of  the  packet  activity.  Such  observation 
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can  be  exercised  by  scanning  a  narrow  field  of  view  (fov)  instru¬ 
ment,  or  using  a  wide  fov  instrument  with  high  resolution  such  as 
an  infrared  imager  [5].  IGWs  with  period  of  15  minutes  to 
several  hours  have  horizontal  wave  length  of  30  to  100  km,  which 
translates  into  the  span  half  angle  of  10  degrees  (30  km)  to  30 
degrees  (100  km)  for  vertical  observation.  At  least  two  data  are 
needed  to  detect  a  wave;  therefore,  the  resolution  needs  to  be  at 
least  half  of  the  span  angle.  Because  there  must  be  at  least  two 
spatial  samples,  the  entire  observation  needs  to  occur  at  more 
than  10  to  30  degrees  half  angle.  This  translates  into  the 
recommended  use  of:  a)  scanning  an  instrument  with  the  fov  5 
degrees  half  angle  or  less  scanning  30  degrees  or  larger  half 
angle,  b)  use  of  an  instrument  with  30  degrees  half  angle  or 
larger  fov,  with  a  spatial  resolution  of  5  degrees  half  angle  or 
less.  The  numbers  for  fov  need  to  be  smaller  for  off -zenith 
observations  according  to  the  factor  cosfl0 ,  where  80  is  the  angle 
of  the  line  of  sight  with  respect  to  the  zenith. 

5.  A  more  theoretical  study  of  the  maximum  entropy  method  is  needed 
to  overcome  some  of  the  problems.  Its  capability  of  working  with 
very  small  amounts  of  data  without  the  need  for  a  window  makes  MEM 
suitable  for  the  signal  analysis  of  airglow  oscillation  data. 
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The  Internal  Gravity  Vaves 

In  this  section,  the  theory  of  the  atmospheric  internal  acoustic/ 
gravity  wave  is  described  based  upon  a  general  theory  of  waves  in 
fluids  [13],  [18],  In  an  atmospheric  gravity  wave,  both  compressional 
and  gravitational  forces  (buoyancy)  cause  wave  propagation.  At  the 
high  frequency  end  of  the  spectrum,  these  waves  approach  to  acoustic 
waves,  while  at  the  low  frequency  end  approach  to  pure  gravity  waves 
such  as  the  water  surface  wave. 

7he  basic  properties' of’ tlie.  acoustic  wave  are:  1)  Isotropic  2) 

* 

Non-dispersive.  The  speed  of  sound  in  the  air  "c,"  for  example,  is 
given  by: 


c  -  J  y?  /  p  -  J  7RT 


where 


(A.l) 


P  -  atmospheric  pressure 
p  -  atmospheric  density 
R  -  gas  constant 
7  -  Cp/Cv  a  1.4 


(N/m2) 

(kg/m3) 

8.314  (J/mol.K) 

(for  air  below  100  km) 


There  are  two  types  of  propagating  gravity  waves:  1)  surface 
waves  and  2)  internal  waves.  An  example  of  a  surface  wave  is  a  wave 
on  the  ocean  surface.  This  type  of  wave  can  exist  only  at  a  boundary 
or  discontinuity  in  the  density  of  the  medium.  This  wave  is  isotropic 
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in  Che  horizontal  direction  which  is  the  propagating  direction,  since 
the  effect  of  gravity  is  Che  same  for  all  directions.  However,  unlike 
acoustic  waves,  surface  waves  are  dispersive.  For  example,  the  wave 
velocity  of  a  long  period  sea  wave  "cs"  is  given  by: 

cs  -  J  g.\x  /  2 ir  (A. 2) 

where  g  is  the  acceleration  o:  gravity  and  Ax  is  the  horizontal  wave¬ 
length. 

Internal  gravity  waves  can  exist  in  the  atmosphere  where  the 
density  decreases  continuously  with  increasing  altitude.  This  type  c: 
wave  can  freely  propagate  to  all  elevation  angles,  and  is  not  limited 
to  horizontal  propagation  as  are  surface  waves.  Since  the  restoring 
force  due  to  gravity  is  always  in  a  fixed  direction,  there  is  no 
reason  for  isotropy  of  propagation.  Hence  internal  gravity  waves  are 
anisotropic  and  dispersive. 

Properties  of  Internal  Acoust ic- Gravity  Waves 

The  simplest  model  is  based  upon  the  assumption  that  the  atmos¬ 
phere  is: 

a)  continuously  stratified  and  compressible, 

b)  isothermal  and  uniform  in  composition, 

c)  stationary  in  the  absence  of  waves, 

d)  invisid  (no  dissipation  of  energy  in  the  gas). 

Because  the  actual  atmosphere  has  non-uniform  and  variable  nature, 
such  model  is  not  always  appropriate. 
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This  model  also  assumes  chat  the  wave  motion  is  adiabatic  and 
produces  small  variations  with  a  period  of  less  than  one  day.  Exact 
solutions  are  found  for  two  dimensional  plane  waves  of  the  form 
exp(i(wc-kxx-kzz) ) ,  where  u  is  the  angular  frequency,  kx,  k2  are  the 
real  horizontal  and  vertical  wave  numbers.  This  form  of  waves  can  be 
the  solution  when  these  parameters  satisfy  the  following  dispersion 
equation  [ 13 ] : 


-  w^c^(kx+k|)  +  <7-l)g^k|  - 


2  2  2 
7  gw 

4c^ 


-  0 


(A.  3) 


This  can  be  rewritten  as: 


-  w2(k2  +  k 1)  +  w|k2 


-  0 


(A. 4) 


where 


wA  “  7g/2c 


wg  “  Jl  -  1  g/c 


(A. 5) 


The  properties  of  the  wave  motion  can  be  deduced  from  this 
relationship.  For  any  pair  of  wave  numbers  (kx,kz),  two  distinct 
values  of  w  exist.  Their  asymptotic  roots  are: 


^2  -  c2(kx  +  k2) 


(A. 6) 
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k|+k2 


(A.7) 


Equation  (A. 6)  represents  the  ordinary  relationship  of  the  non* 
dispersive  property  of  sound  waves,  and  is  valid  when  The 

equation  (A.7)  represents  the  dispersion  relation  of  an  internal 
gravity  wave  under  the  asymptotic  condition  ut  «  u%. 

The  wave  can  propagate  only  if  kx  and  kz  are  positive,  which  is 
true  when  u  <  u9  .  Hence,  internal  atmospheric  gravity  waves  can  exist 
only  at  periods  greater  than  a  certain  value  Tg  which  is  defined  as: 

T,  -  2irc/g/7~l  (A.  8) 

The  value  Tg  is  known  as  the  Brunt *Vaisala  period.  This  is  the 
natural  resonance  period  of  a  displaced  parcel  of  air  under  the 
restoring  force  due  to  buoyancy  alone. 

Similarly,  internal  acoustic  waves  require  w  >  uA .  Thus,  internal 
atmospheric  acoustic  waves  can  exist  only  at  periods  smaller  than  a 
characteristic  value  T^>  which  is  defined  as: 

T^  -  Airc/yg  (A.  9) 

The  value  TA  is  known  as  the  acoustic  resonance  period.  This  is  the 

f  • 

resonant  period  of  the  whole  atmosphere  with  the  restoring  force  due 
only  to  compression. 
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Since  the  ratio  7  -  Cp/Cv  is  less  t  an  2  within  the  homosphere 
then  <  Tg  and  the  two  types  of  internal  waves  are  completely 
decoupled.  Fig.  A-l  shows  the  permissible  values  of  k*  and  kz  for 
several  different  periods.  The  acoustic  waves  form  a  family  of 
ellipses  and  the  gravity  wa\u  form  a  set  of  hyperbolas.  At  periods 
less  than  one  minutes,  the  ellipses  are  almost  circular  implying 
isotropic  wave  motion,  or  the  regular  acoustic  wave,  i.e.  sound. 

Internal  gravity  waves  usually  travel  much  slower  than  the  speed 
of  sound.  In  general,  the  shorter  the  period,  the  smaller  the  phase 
speed  and  the  shorter  the  wavelength  of  the  wave.  Internal  gravity 
waves  also  have  a  property  that  the  direction  of  energy  flow  is 
different  from  that  of  the  phase  propagation.  For  a  long  period 
wave,  both  directions  are  almost  normal  to  each  other. 

For  shorter  period  gravity  waves,  the  motion  of  the  air  parcel  is 
not  only  transverse  but  is  slightly  elliptical,  resulting  from  the 
combination  of  the  longitudinal  compressional  forces  and  transverse 

4 

buoyancy  forces. 

As  a  consequence  of  the  conservation  of  energy  flux  in  an  invisid 
atmosphere  where  density  decreases  exponentially,  amplitude  of  waves 
increases  exponentially  with  altitude. 

Fig.  A- 2  shows  the  relationship  between  the  horizontal  wavelength 
and  phase  velocity  for  several  phase  propagation  angle  theta  (0)  of 
IGW.  The  direction  of  phase  velocity  is  given  by  9  and  shown  to  be 
downward.  Using  the  equation  (A. 7),  6  satisfies: 
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kj(m-l) 


The  family  of  ellipsoids  represents  the  acoustic-wave  solutions  to  the 
equatidri  ( A . 6 )  The  family  of  hyperbolas  represents  the  IGMs, 
solutions  to  equation  (A. 7).  Each  curve  is  marked  by  oscillation 
periods  in  minutes. 


Fig.  A-l.  Permissible  IGWs 
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Phase  Velocity  (m/s)  Phase  Velocity  (m/j)  Phase  Velocity  (m/4) 
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Horizontal  Wavelength  (km) 


Fig.  A- 2.  1GW  topography 


195 


(w2  -  u|) 

fG  -  “  w2  -  u2(k2  +  k|)  +  u|k2 

c2 


is  the  left-hand  side  of  the  dispersion  equation  (A. 4). 
(A. 13)  into  (A. 11)  and  (A. 12)  yields: 

wkx  (w2  -u2 ) 


a>4 

7  - « 

c2 


w3kz 


W4 

7  -  44 

c2 


This  is  the  relationship  between  the  packet  velocity 
wave  parameters.  The  elevation  angle  of  the  wave  energy 
q,  may  be  specified  with  the  following  equation: 


kz  w2  Pz 

tan  a  - - -  — 

kx  "2  -  4  P* 


(A. 13) 

Substituting 


(A.  1-0 


( A . It) 

(18j  and  the 
propagation , 


(A. 16) 
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aendix  B:  Rotational 


aerature 


After  the  optical  emission  spectrum  is  obtained,  the  band  emission 
intensity  and  the  rotational  temperature  for  a  specified  band  are  to 
be  calculated.  First,  the  definition  of  the  rotational  temperature  is 
presented. 

The  hydroxyl  can  be  modeled  as  a  dipole  consisting  of  two  mass 
points  with  positive  and  negative  charge  connected  by  a  "massless 
spring"  as  in  Fig.  B-l.  When  generated  by  the  chemical  reaction: 


03  +  H  - >  OH*  +  02 


(B.l) 


where  i  denotes  the  vibrational  level.  Energy  for  hydroxyl  in  this 
reaction  is  3.3  eV  [10]  which  is  too  small  for  the  electronic  state  to 
be  excited,  so  they  remain  at  the  ground  state.  This  energy  is  large 
enough  to  excite  the  vibrational  level  whose  quantum  number  is  less 
than  or  equal  to  nine.  These  vibrational  levels  are  associated  with 
rotational  levels .  These  rotational  levels  cause  the  rotational 
emission  lines  in  the  emissive  vibrational  level  transitions.  As  a 
result,  those  transition  emissions  produce  emission  bands  that  are 
made  by  numerous  rotational  lines.  The  distribution  of  emission 
intensity  among  the  rotational  lines  within  a  single  band  is  a 
function  of  temperature.  In  this  context,  "temperature"  is  the 
kinetic  temperature  of  hydroxyl  that  can  be  assumed  identical  to  the 
kinetic  temperature  of  the  major  species  as  long. as  the  time  scale  of 
interest  is  larger  than  the  inverse  of  the  collision  frequency  which 
is  on  the  order  of  few  milliseconds.  With  this  method,  the  tempera- 
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ROTATION 


Fig.  B-l.  Hydroxyl  model  (lCf  . 


ture  can  be  calculated  from  the  distribution  of  several  ro  esc  ion  Ci  1 
line  intensities.  Theoretically,  only  two  rotational  lines  are 
necessary  and  sufficient  to  calculate  the  temperature.  But  with  veal 
data  this  method  is  too  sensitive  to  any  error  in  the  line  intensities 
including  the  effect  of  noise.  The  theoretical  ratio  of  rotational 
lines  in  a  band  for  given  hydroxyl  concentration  and  temperature  can 
be  calculated  as  follows.  Let  v'  and  v' '  denote  the  initial  and  final 
vibrational  states,  thus  making  (v',v* ')  band.  Let  j'  and  j''  denote 
the  initial  and  final  rotational  states.  The  concentration  of  the 
initial  state  is  given  as  Nv'  .  Then  the  concentration  of  each 
rotational  state  j',  Nv'j’  is  represented  as  [9]: 
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Nv'j ' 


^  (2j'+l) 
Qr 


exp(- 


Bi'd'+Dhc 

kT 


(B.2) 


B  -  — ■£—  ~  16.28  (cm" 1 )  [9] 

8jt2  cl 


(3.3) 


Qr 


03 


S  (2j+l) 
i*o 


,  Bi(i+l)hc, 

exp( - - ) 

kT 


(B.4) 


where  T  is  the  temperature,  Qr  is  the  partition  function,  B  is  the 
rotational  constant,  h  is  the  Planck  constant,  c  is  the  speed  of 
light  and  k  is  the  Boltzmann  constant.  The  emissive  transition 
occurs  only  according  to  the  selection  rule.  In  this  case  the 
selection  rule  is  Aj  -  O.+l;  the  change  of  the  rotational  level  is 
either  none,  plus  one  or  minus  one,  corresponding  to  what  are 
customary  called  Q,  P  and  R  branches.  Each  branch  is  divided  into 
two  smaller  subsets..  This  is  based  on  the  two  different  total 
electronic  angular  momentum  about  the  internuclear  axis.  It  is 
represented  by  a  quantum  number  Cl,  that  can  take  either  1/2  or  3/2 
[10].  Each  branch  of  P,  Q  and  R  lines  is  denoted  with  a  subscript  of 
1  or  2  depending  whether  H-3/2  (Pi,  Ql  and  Rl )  or  1/2  (P2  ,  Qz  and  Rz ) . 
The  most  straightforward  method  to  calculate  the  rotational  tempera¬ 
ture  is  to  use  the  rotational  lines  of  a  single  branch.  In  case  of 
(3,1)  band  of  hydroxyl,  the  P  branch  is  commonly  used.  This  is 
because  the  Q  branch  lines  are  too  closely  spaced  for  the  interfero¬ 
meter  spectrometer  to  resolve,  and  the  R  branch  falls  into  the 
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wavelengths  chat  are  in  the  water  vapor  absorption  region.  in  prac¬ 
tice,  more  than  three  lines  are  used  to  calculate  the  rotational 
temperature.  Among  those  are  the  methods  developed  by  Hill,  Espy  and 
Niple  independently  [31],  [35]  using  the  least  square  fitting  of  a 
synthesized  spectral  lines  of  one  band  to  the  experimental  data  of  die 
spectral  lines  in  the  same  band. 
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Appendix  C:  Linearization  of  the 
Adiabatic  State  Equation 

This  derivation  of  the  linearized  equation  of  adiabatic  change  of 
state  is  adopted  from  Hines  {11] .  This  gives  an  important  expression 
for  the  derivative  of  the  vertical  velocity  component  w  of  v  - 
[u,0,w]  with  the  assumed  form: 

w  -  W(z)exp(i(wt-kxx))  (C.l) 

The  equation  of  horizontal  acceleration  yields: 

iwp0u  —  ikxp'  (C.2) 


or 


p'  7«u 

po  kxCz 


(C.5) 


with  the  speed  of  sound  C-(7po/po )**,  po  and  po  being  the  unperturbed 
mass  density  and  pressure  and  those  with  apostrophe  being  their 
perturbations.  The  Lagrangian  derivative  of  the  pressure  is: 


I 


Dp  3p '  dpo 

—  -  —  +  w  — 

Dt  3t  dz 


iup' 


Po 
w  — 
H 


(C.4) 
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•  1  , 

where  H-Hp-  -(d(ln  po)/dz)  -  po/pog  -  Cz/7g  is,  the  pressure  scale 
height,  g  being  the  gravitational  acceleration.-  The  equation  of 
adiabatic  change  of  state  can  be  represented  as: 

7  Dp  1  Dp 

~7  Dc  ”  ~p~  DC  (C'3) 

rO  r  0 


and  the  continuity  equation: 

1  Dp  ow 

-  —  -  -Vv  -  ikxu  -  — 

P0  Dc  ~  3 z 


(C.6) 


With  (C . 5)  and  (C.6),  (C.3)  yields: 


ikxu  - 


dw 

dz 


1  3w 
7H  3z 


and  so: 


ikxu 


k*C2 


k2C2  -  W2 


3w  _w  , 
dz  ~  7H  ' 


From  (C.6)  and  (C.8), 


Vv 


-  ikxu  + 


dw 

dz 


~  ^C2 
k2C2-*2 


.  dv  w 
'  dz  ~  7H  ' 


(C.  7) 


(C.8) 


201 


202 


-  - * -  (  -  kjgw  >  (C.9) 

*lc  82 

The  adiabatic  equation  of  state  can  also  be  written  in  the  following 
form: 

22  -  iui*  +  w  ^22.  .  .(.7-111°  E£  (C.10) 

Dc  dz  po  Dt 

(C.6)  and  (C.10)  yield: 

To  (z)  +  T’(z)  -  To(f)  (1  +  i  ^  Vv  )  (C.ll) 

w  If 

This  is  the  equation  (2.73).  The  associated  equation  for  the  change 
in  density,  (2.72),  is  derived  as  follows.  From  the  continuity 
equation,  the  following  form  is  derived: 

iwN '  -  -  No  Vv  -  iwh  ^  (C.12) 

oz 


or 


N'  -  -  Na  Vv  -  h  dS2 

iw  ~  dz 


(C. 13) 


since  VN  -  dNo/dz,  or  No  is  constant  for  the  same  altitude  with 
different  x  or  y  coordinates.  Equation  (C.13)  leads  to: 
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No(z)  + 


''(*>  ‘  No<f)+hi5fn  +  '-S^-hf!!)|C 


-  No(0  +  h  +  t  »°  *»)f-  h  ^ 


No(f)(  1  +  A  Vv  >._ 

w  |  f 


(C.U) 


which  is  che  equation  (2.73). 


